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Abstract 

A propeller which operates in the shear flow found near the aft end of marine vehicles experiences 
an intimate coupling between the propeller’s induced velocity field and the rotational inflow. The 
presence of the propeller's induced velocity field causes the inflow to accelerate, which redistributes 
the vorticity present in the inflow. This redistribution causes a change in the nominal propeller 
inflow. Because the propeller now experiences a different inflow, the propeller induced velocity field 
is altered. Thus, there is an intimate coupling between the vorticity present in the fluid inflow and 
the propeller-generated induced velocity field. 

Lifting surface propeller blade design codes are incapable of analytically representing the vortical 
interaction between the induced velocity field and the rotational inflow. An additional difficulty is 
encountered as downstream blade elements, which would be present in a multi-component propulsor, 
pass through the singularity wake sheets shed by upstream components. For these reasons, the 
propeller blade design code should be coupled with an external flow solver which is capable of 
transporting vorticity. 

Previous researchers have coupled propeller blade design codes with Reynolds Averaged Navier 
Stokes (RANS) flow solvers. This powerful method made possible multi-element propulsor design in 
the presence of rotational inflow. Unfortunately, the use of a RANS code is costly in terms of time 
and computational resources. 

This thesis focuses on coupling a propeller blade design code with an axisymmetric, multi- 
element through-flow code, developed by Drela. The throughflow code uses an integral boundary 
layer method to solve for the boundary layer flow, and a streamline curvature formulation to solve 
for the in viscid, outer flow. The main advantage of the present method over previous methods is an 
order of magnitude reduction in computation time. Validation cases were performed to validate the 
various components of the coupling procedure, as well as the coupling methodology as a whole. A 
design case is presented which shows the use of this methodology in design. 

Thesis Supervisor: Justin E. Kerwin 
Title: Professor of Naval Architecture 
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Chapter 1 



This Thesis in the Propeller 
Design Process 



A propeller which operates in the shear flow found near the aft end of marine vehicles experiences 
an intimate coupling between the propeller’s induced velocity field and the rotational inflow. The 
presence of the propeller’s induced velocity field causes the inflow to accelerate, which redistributes 
the vorticitv present in the inflow [14]. This redistribution causes a change in the total inflow 
velocity to the propeller. Because the propeller now experiences a different inflow, the propeller 
induced velocity field is altered. Thus, there is an intimate coupling between the vorticity present 
in the fluid inflow and the propeller-generated induced velocity field. To design more efficient, and 
more complex, propulsor geometries, it is necessary to accurately account for this physical interaction 
phenomena. 

A potential, vortex lattice lifting surface propeller blade design code such as Kerwin’s Propeller 
Blade Design (PBD) code is incapable of analytically representing the vortical interaction between 
the induced velocity field and the rotational inflow. An additional difficulty is encountered as 
downstream blade elements, which would be present in a multi-component propulsor, pass through 
the singularity wake sheets shed by upstream components. For these reasons, the propeller blade 
design code should be coupled with a throughflow fluid solver which is capable of modelling vorticity 
transport. 

This thesis focuses on coupling a propeller blade design code developed at MIT by Kerwin [5] 
with an axisymmetric, multi-element throughflow code, developed at MIT by Drela [3]. Drela's 
throughflow solver uses an integral boundary layer method to solve for the boundary layer flow, 
and a streamline curvature formulation to solve for the in viscid, outer flow. The main advantage of 
the streamline curvature throughflow code is an order of magnitude reduction in computation time 
compared with the RANS coupling methodology [14]. 



1.1 Fluid Velocity Terminology 



Before proceeding further, it is necassary to first clear up what the propeller designer means when 
talking about fluid velocities. The nominal inflow is the fluid velocity in the region of the propeller 
with the propeller not operating. The total inflow is the fluid velocity with the propeller operating. 
The total inflow velocity is composed of an induced velocity, which is induced by the presence of the 
propeller blade and trailing wake singularity distributions, and an effective inflow. 



^ total — 



f induced T ^ effective 



( 1 - 1 ) 



Figure 1-1 graphically shows the relative magnitudes and directions of the total, induced, and 
effective inflow velocities for a typical propeller. 

While the nominal inflow is the inflow field when there is no propeller present, it is not equal 
to the effective inflow. This is due to the aforementioned interaction effects between the propeller 
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Total Inflow Induced Velocities 




Effective Inflow 




Figure 1-1: The total, induced, and effective inflow velocity fields generated by Propeller 4119, 
J=0.833, in inviscid flow. 

vortex lattice and the vorticity in the inflow. This can be seen from the definition of the vorticity 
vector. 



J = ( 1 . 2 ) 

A non-radially-uniform meridional inflow suggests the presence of a tangential component of vortic- 
ity. 



dV r dV x 

dx dr 



(1.3) 



The tangential vorticity can be thought of as a ring which circumscribes the propeller shaft [18]. 
The propeller induced velocity causes a contraction of the tangential vorticity ring. As this ring 
contracts, Kelvin’s thereom requires that the strength of the vorticity remain constant. Therefore, 
from equation (1.3), it is evident there is a change in to counteract the change in If there 
were no tangential vorticity present in the inflow, though, the nominal and effective inflows are 
exactly equal. This fact is used later to validate the coupling methodology in this thesis. 



1.2 Effective Inflow 



All propeller design and analysis takes place in the effective inflow. Recall that the effective inflow 
is the total inflow with the effects of the propeller (*.e., the propeller induced velocity component) 
removed. The effective inflow is a non-physical flow field. It can not be measured in a water tunnel 
or behind a ship in operation. And, yet, determining the effective velocity field is critical to the 
success of the propeller design process. 

Historically, to determine the effective inflow velocity, a model of the ship is constructed, and 
velocity profile surveys are taken in the propeller disk plane of the unappended model. This is the 
nominal inflow. The effective inflow velocity is taken as a fraction of the nominal inflow. 



Veff = 



n -w T \ 



(1.4) 



As shown above, though, this linear scaling across all radii means that in the hub region, where 
the boundary layer dominates the inviscid flow, the velocities are too high, and the velocities over 
the rest of the blade span are too low. In fact, this method produces the correct velocity at only 
one radius along the propeller span! 

A more modern method is to computationally compute the flow in the region of the propeller by 
solving the equations of fluid motion, accounting for the presence of the propeller. Thus, the exact 
total inflow velocity is known at every radii across the propeller span. Once the induced velocity is 
subtracted from the total inflow velocity, the effective inflow is known. 

Besides the effective inflow velocity, the propulsor designer also requires an accurate knowledge 
of the thrust power which must be generated by the propulsor at the design speed. 
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AY, A ’q,t] 

Propeller Diameter 
Spanwise Circulation, T(r) 
Propeller RPM 
Cavitation Index cr 



Table 1.1: Outputs from a lifting line analysis. 



1.3 The Coupled Hull Flow Resistance Problem 

To accurately design a propeller requires an accurate knowledge of the ship’s resistance. The in- 
creased fluid velocity along the stern of the vessel due to the presence of the operating propeller 
causes an increase in drag since there is more incomplete pressure recovery. The vortical interaction 
also wreaks havoc upon the traditional tow tank resistance calculation, since, for vessels with high 
aft prismatic coefficients, the presence of the propeller can inhibit boundary layer separation. 

Historically, the values for the resistance of the ship with and without a propeller operating are 
related through the thrust deduction coefficient. 

R T = (l-i)T (1.5) 

Modern computation techniques which model the interaction of the propeller and the fluid flow 
around the hull solve both the effective inflow and thrust deduction problems by nearly exactly 
solving for the fluid flow over the hull with the propeller operating. In this method, the entire 
submarine and duct, if present, is modelled in the computer. The shear flow is exactly computed 
along the body, and the correct propeller interaction effects are also modelled. 

Once the effective velocity field is known, and the ship, or submarine, resistance is known, we 
can get on with the task of designing a propeller. 



1.4 Modern Design Techniques 

The starting point for the propeller design process is the effective inflow. The designer then uses 
a lifting line theory to optimize the gross propeller characteristics. In lifting line theory the blades 
are represented as two dimensional “sticks” of bound vorticity [13]. Refinements upon the basic 
lifting line model allow for hub and duct modelling, and multiple blade rows [2]. Bowles’ advanced 
lifting line method [2] also allows for the cavitation and strength considerations in his lifting line 
development. Not included in this analysis are skew and rake effects. The outputs of the lifting line 
analysis shown in Table 1.1 give the several gross propulsor characteristics. 

Once the lifting line analysis is complete, a lifting surface propeller code is used to find the 
two dimensional blade shapes, pitched about a raked and skewed genarator line, which produce the 
required circulation F(r). While further discussion of lifting-surface theory is delayed until Section 
1.5, there are several key output parameters from lifting surface design. 

In the serial progression from lifting line to lifting surface design shown in Figure 1-2, the effective 
inflow is assumed constant. However, a quick check of the numerically predicted induced velocity 
components against the empirically assumed induced velocity components (remember, the designer 
had to scale a physical total inflow to get to the effective inflow by subtracting the induced velocity 
components) would show that there are differences. This should drive the propeller designer to 
rederive a new effective inflow and optimize a new propeller. Another tact is to couple the preliminary 
design tools with a viscous ( i.e ., able to model the viscosity in the governing fluid equations of 
motion) or inviscid Euler throughflow solver which can accurately model the effects of the propeller 
vorticity upon the global flow solution. The requirement for preliminary design, however, is a fast, 
computationally “cheap” solution to allow proper exploration of the design space. The current 
RANS coupling is not fast in that sense. As will be shown later, an axisymmetric streamtube solver 
meets both requirements. 
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Lifting Line Design Blade Representation 



Lifting Surface Design Blade Representation 
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No current coupling 
between a lifting line 
blade design code and 
throughflow solvers due 
to computational expense 
of RANS-based solvers. 



Current Design 
Techniques Iterate 
A Lifting Surface 
Design Code with 
RANS-based 
throughflow solvers 



Numerical Representation of Body and Propeller Forces in Throughflow Fluid Solver 




Figure 1-2: The prelimary propeller design process relies upon ail assumed effective wake 
( i. e., velocity field) as input to a lifting line analysis to optimize the propeller and lifting surface 
analysis to design the three dimensional blade shapes. There is a limited amount of coupling with 
throughflow solvers to iterate upon the effective wake. 
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Performance 
Blade Offsets 



Spanwise Circulation 



Local Pressure Coefficient 



AY, 10 Aq, t /, to check thrust and torque. 
Expressed in the traditional propeller design pa- 
rameters of rake, skew, pitch, camber, and thick- 
ness distributions this specifies the geometry of the 
blade for construction. 

The final circulation distribution may not match 
the original desired circulation distribution due to 
three dimensional rake and skew effects. 

This gives and indication of cavitation. 



Table 1.2: Outputs from a lifting surface analysis. 



1.5 Lifting- Surface Propeller Blade Design and Analysis 

In a traditional propeller lifting surface propeller code, a grid lattice is placed on the blade mean 
camber surface, the hub and duct (if present) and the trailing wake system [15]. Each lattice segment 
is assigned a strength of vorticity. On the solid body surfaces, such as the blade and hub, a control 
point is placed near the center of the grid lattice. The strength of each vortex lattice segment 
is assigned by satisfying the kinematic boundary condition that the flow must be tangent to the 
propeller, hub, and duct surfaces at every control point. 

Mathematically, the propeller problem involves a simple matrix equation. By attacking the 
geometry of the problem, an influence matrix is formulated which gives the velocity induced by 
a unit strength of vorticity along every vortex lattice upon every control point - an [7A r F]luence 
matrix. When multiplied by the actual vortex segment strengths, the induced velocity at every 
control point is known. 



[INF] 



= Vin 



duced 



(1.6) 



/noindent To model the effects of blade thickness, a souce distribution is placed coincident with 
the blade vortex lattice system. 

Next, to satisfy the kinematic boundary condition, the physics of the problem dictate that the 
component of the induced velocity normal to the blade, when added to the component of the effective 
inflow velocity normal to the blade must be zero. 



[[INF] [f] + Veff 



71=0 



(1.7) 



Equation (1.7) is the heart of the propeller lifting surface code. 

A lifting surface propeller code can be used for either the design of a new propeller or an analysis 
of an exsisting propeller. In the case of a propeller design, a radial loading distribution is prescribed, 
and a blade shape is found which produces the desired loading by manipulating the geometry of 
the blades. In blade analysis, the unique strength of the vortex lattices segments is found which 
satisfies equation (1.7). From equation (1.7) it is seen that an incorrect effective inflow velocity will 
give either an incorrect blade shape, or an incorrect prediction of blade loading. 

The propeller forces resulting from the vorticity and source distributions are calculated from 
the Kutta-Joukowski and Lagally theorems, respectively. A Lighthill leading edge suction force 
correction is applied to these forces, and the propeller’s sectional drag is calculated either based on 
strip wise two dimensional empirical drag coefficients, or a stripwise 2-D integral boundary layer 
calculation. 

The propeller blade analysis code used in this study is an extension of a previously reported 
lifting surface design code [12, 5]. The extensions replace the image hub and duct with vortex lattice 
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Propeller 41 19 Vortex Lattice Representation 




Figure 1-3: 3 bladed propeller 4119 oil a straight shaft. Notice that the hub discretization extends 
upstream of the hub. The inner radius of the transition wake vortex lattice is the innermost hub 
streamline. 



representations of the hub and duct [14, 17]. The computer representation of the blade and hub 
vortex lattice, and transition wake is shown in Figure 1.5. 

The end result of the propeller blade analysis routine is the complete circumferential mean 
induced velocity over the blade surface. The circulation, T# ,is computed from the circumferential 
mean tangential induced velocity by the following relationship which follows from the Euler turbine 
equation. 



ffl(r) 



-2kvV; 

z 



( 1 . 8 ) 



Equation (1.8) is at the heart of the present coupling technique with the streamline curvature 
method. In the streamline curvature method, as in turbomachinery mean line design, the designer 
is most interested in determining the amount of turning, or swirl, being applied to the flow through 
the action of the blades. Therefore, in the coupling with the streamline curvature method, the swirl 
7*(A\ 7 0) is directly calculated by multiplying V 0 * by the local radius r. The use of the swirl velocity 
to couple between the propeller code and the throughflow solver is unique to the present research. 
In previous coupling methodologies, the propeller is represented as distributed body forces within 
the RANS throughflow solver. 



1.6 Coupling with a Viscous Throughflow Solver 

The major shortcoming of the lifting surface formulation described above is the inability to model 
the vorticity interaction. As previsouly mentioned, the propeller induction velocities cause a redis- 
tribution of the vorticity present in the inflow. This interaction is critical in the design of modern, 
full stern submarines, where the stern prismatic coefficient is so high that the propulsor actually 
inhibits separation of the boundary layer. 

The second area of propulsor design affected by this shortcoming is in the design of multi-element 
propulsors. The problem is that at points where the trailing wake vorticity from the upstream blade 
rows intersect the control points of the downstream blades, the mathematics of the problem produces 
a singular solution. Multi-element propulsor design is important in the design of modern submarine 
propulsors, and in the design of waterjet pumps. 

And yet, the lifting surface formulation is unequalled by any other technique in its ability to 
rapidly and accurately predict the forces produced by a propeller. But it can not model vorticity 
interaction. Modern throughflow solvers, though, can accurately model vorticity interaction, but are 
too computationally expensive to model the whole propeller problem in three dimensions. Modelling 
the three dimensional propeller blade problem in a fast vortex-lattice lifting-surface code means 
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that a computationaly axisymmetric “cheap” viscous throughfiow solver can be used to solve for the 
circumferential mean flow over the body. 

Past researchers have coupled propeller design codes with axisymmetric RANS throughfiow 
solvers. This method is extremely robust, but also extremely time and computationally inten- 
sive. Thus, it seems logical to couple the propeller blade design code with a streamline curvature 
code, because the streamline curvature code runs extremely quickly. A simple test case, such as a 
three bladed propeller on a straight shaft in inviscid flow, takes around 12 hours to run on a typical 
work station. More complex cases, such as a submarine propulsor, take on the order of days. As a 
comparison, it takes less than 15 minutes to run the simple propeller in inviscid flow with the present 
streamline curvature coupling method. Obviously long computer run times do not leave much tune 
for the engineer to refine and optimize the design since so much time is taken up by the analysis of a 
single design point. The potential time reduction would seem well worth the effort of incorporating 
the streamline curvature method into the preliminary design as it would allow for more exploratory 
design efforts, rather than computational analysis. 
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Chapter 2 

Flow Theory 



2.1 Computational Solution to Fluid Flow Problems 

The governing equations of motion for any fluid flow are the Navier-Stokes Equation [22], written 
here in compact vector form. 



^ = --VP+ t/V 2 r + VH (2.1) 

Dt p 

The solution to Equation (2.1) requires discretization in time and space. To ensure accuracy and 
stability of the numeric solution, the discretization sizes are directly related to the smallest length 
and time scales present in the flow. It can be shown that due to the extremely small length and time 
scales involved in boundary layer flows, a complete, numeric solution to the Navier-Stokes equations 
would require computational resources far in excess of those available today. Therefore, various 
levels of approximations are introduced to find solutions to equation (2.1). 

To make the Navier-Stokes equation tractable, one fundamental approximation is to ignore the 
viscous terms. This results in the Euler equations of motion for an inviscid fluid. Hydrodvnamicists 
are able to simplify the problem further by considering the fluid as incompressible. For extremely 
simple flow cases, there exist numerous, classical theoretical solutions to the incompressible, inviscid 
flow equations. The value of these classical results is not to be underestimated in that they provide 
an irrefutable answer by which numeric solution techniques can be validated. 

Solutions to real world flow problems are only possible using numerical schemes. The large 
success of incompressible, inviscid numeric computer codes shows that even with the inviscid and 
incompressible approximations, the majority of the physics of the flow is still correctly modelled. 
There are, however, ways to introduce the effects of the viscosity without the computational expense 
of solving the fully three dimensional Navier-Stokes equations. 

Since viscous effects are mostly confined to a small region of the flow domain, an effective ap- 
proximation is to use a “divide-and-conquer’’ strategy. In this technique, an inviscid Euler solver 
is used to solve the majority of the flow domain, and a boundary layer solver is used to solve for 
the boundary layer growth along a solid wall. The two codes co-exist by matching the boundary 
conditions (pressure and velocity) at the interface between the two flow regions. This type of cou- 
pling is a displacment body type of scheme and is the scheme employed by the numeric codes used 
in this thesis. The basic idea is to use a fast, inviscid solver to solve for the inviscid flow and a fast 
boundary layer routine to solve for the the boundary layer flow. 

2.1.1 The Boundary Layer Equations 

A significant approximation to the full Navier Stokes equation can be made through the application 
of dimensional analysis [25]. The resulting simplfied equations, usually termed the Thin Shear Layer 
Equations or the Boundary Layer equations are quite simplified indeed. They are written here with 
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the Conservation of Mass equation for completeness in the subsequent derivation of the Integral 
Boundary Laver Equations. 



Continuity §7 + §£ 

Momentum + t'lp- 

ar dx ay 



= 0 








( 2 . 2 ) 



(2.3) 



To formulate the so-called integral boundary layer equations. Equation ( 2 . 2 ) is pre-multiplied 
by the quantity (u — u e ) and subtracted from Equation (2.3) [25]. The resulting equation is then 
integrated over the domain for the boundary to infinity. 



t 0 0 f°° du 

-j- = — J (u e -u)dy+ — J I u(u e - u)dy+ J (u e - u)dy - u e v w (2.4) 

where r shear stress on the solid boundary 

u t inviscid edge velocity 

u(y) local velocity 

v w velocity through the wall, usually termed as wall 

transpiration 

The key at this point, is to recognize the definition of the displacement thickness and momentum 
thickness buried beneath some algebra within Equation (2.4). Recall the following equations are the 
definitions for the displacement and momentum thicknesses. 



Displacement Thickness S* 

Momentum Thickness 0 

Inserting Equations (2.5) and (2.6) into Equation 
boundary layer equations. 



u 

(1 )dy 

0 u e 


(2.5) 


u u 

_(1 -~)dy 

0 u t u e 


( 2 . 6 ) 



(2.4) gives the common form of the integral 



T 



pin 




1 d dO 1 du e 

) + ~j — l-(20-f(5 ) — 

u - at dx u e dx 



u e 



(2.7) 



Ignoring the wall transpiration term, and allowing only steady flows, Equation (2.7) reduces still 
further. 



£' = ^ + (2+w) i^ 

2 dx u e dx 



(2.8) 



where H — V’ niomen t uni shape factor 

From a computation standpoint, Equations (2.7) and (2.8) are superior because they are parabolic 
in nature and can be solved through a marching routine. That is, starting from a given upstream 
starting condition, the solution downstream is found by marching downstream. Equation (2.3), on 
the other hand, is hyperbolic and must be solved globally and, usually, iteratively. 

There remains, however, the complication of a closure relation. In equation (2.8) there are 
more unknowns than equations. Therefore, we must use an auxiliary equation which relates one, 
or more, of the unknowns to the others. This auxiliary equation is termed a closure relation when 
applied to boundary layer flows, and is usually based upon a mix of empiricism and physics. The 
closure relations used within MTFLOW and the PBD-MTFLOW coupling codes are those derived 
by Swafford [23]. 
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2.2 Streamline Curvature Method 



The streamline curvature method is a powerful method for solving for the inviscid flowfield within 
a channel or annular passage [1, 22]. Because this method is inviscid, it must be coupled with a 
boundary layer routine to solve for the boundary layer flow. The two then must operate in tandem 
since the boundary layer tends to displace the inviscid streamlines ofT the body by the displacement 
thickness, which causes a redistribution of the inviscid streamlines. 

The throughflow computer code developed by Drela [3] is based on the conservative formulation 
of the steady state Euler eqations developed by Drela [4]. This approach is an improvement upon 
the original Streamline Curvature Method formulation [19]. Section 2.2.1 shows the basics of the 
streamline curvature method, section 2.2.2 outlines Drela's reformulation of the problem in a finite 
volume sense, and section 2.2.3 shows the MTFLOW computer code implementation, which makes 
another intelligent simplification. 

The streamline curvature method has been previouslv applied to hydrodynamic open propulsors 

[ 20 ]. 

2.2.1 The Original Streamline Curvature Method 

The streamline curvature method is applicable for axisymmetric meridional flow problems. It was 
first derived to solve internal flow problems. In the original streamline curvature method, the 
steady-state conservation of momentum equation is formulated which accounts for changes in radial 
momentum due to meridional curvature, fluid acceleration, and tangential velocity, or swirl, effects 

[ii]- 

Fluid accelerations in the radial direction are due to 

1. Acceleration in the Streamtube Direction. 



Vm 



dVm x 

1>7 J 



which has an outward directed radial component 



/ dv m x . , 

— S1I10 
OS 

2. Meridional Curvature of the Streamtube. The centripetal acceleration is 

tin 

r c 



3. Tangential Velocity. A tangential velocity component (f.e., swirl) implies a centripetal accel- 
eration 



n 



V 



The streamline curvature method is based upon the solution of the equations of motion along a 
streamline in steady flow. In a streamline and normal coordinate system ( s , n) the equations of 
motion are given as the sum of the three acceleration effects cited above[l]. 



dP , T , dV m drV e 

__ + ^=pV m — -pu,— 



(2.9) 



dP , v , 

-- 7 — + P9n = pV K 



( 2 . 10 ) 
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Figure 2-1: Finite Volume geometry and definitions for conservation equations. 



Notice that the integral of equation (2.9) along the streamline reduces to the usual Bernoulli equation 
when the rate of rotation, u>, is zero. 

Equation (2.10) is usually called the radial equilibrium equation. The position of the streamtube 
boundaries directly gives the fluid radius of curvature term k. This simple ODE is integrated in 
the radial direction to derive the mean velocity within the streamtube. The resulting error in the 
conservation of mass equation is used to drive a relaxation procedure to alter the streamlines towards 
their correct positions. 

The main advantage of the streamline curvature method over traditional time-marching methods 
is that itr is orders of magnitude faster, and hence, applicable to industrial design problems [4]. The 
main disadvantage of the streamline curvature routine is that the iterative solvers lack numerical 
stability in regions of supersonic flow, and do not correctly position shocks within the flow. The 
method of Drela [4] was originally developed to correct this deficiency and correctly treat shocks. 

2.2.2 Finite Volume Formulation 

Drela's finite volume throughflow method is based on the streamline curvature method. The steady 
state equations, though, are formulated in a finite volume sense. 



The first governing equation is the conservation of mass. Across any two adjoining finite volumes 



p\VmiSi * ^1 = p?v m2 s 2 • A 2 (2.11) 

Note that in Equation (2.11) that the s\,s 2 vectors lie along the streamline. Conservation of 
momentum is written as the steady state Euler equation ( Navier Stokes equation sans any viscous 
terms ). 
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( 2 . 12 ) 



-Pi j4i 4 (Pi * A\ )i , ml S 1 + II B — P 2^4 2 4 (P2 v m2 s 2 ' ^2) v m2^2 4 FI + B + 



Conservation of energy is formulated in terms of enthalpy. 



7 P\ 
7 “ 1 Pi 
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+ 2* 



2 
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2 

m2 



(2.13) 



All additional constraint imposed is that the average pressures on the faces normal to and along 
the streamtube be equal. 



Pi 4 P 2 = n + + rr (2.14) 

For the subsonic flow cases which are of interest in marine propeller design, the differential equa- 
tions form a set of coupled, elliptic, PDEs. The elliptic nature of Equation (2.11) through Equation 
(2.14) requires an iterative relaxation solution procedure, vice a marching solution technique, which 
would be used in the case of supersonic flow, where the equations are hyperbolic in nature. 

The solution is started by specifying an initial gridline distribution. Each gridline pair is treated 
as a streamtube. Within each streamtube equations (2.11) through (2.14) are solved simultaneously 
at each streamwise stations using a Newton-Rhapson technique. The solution produces the pressures 
FI + and 11“ on the streamtube walls. By calculating influence coefficients for the effects of streamline 
curvature and streamtube area on Afl, a relaxation equation can be used to update the streamtube 
positions. 

2.2.3 Drela’s Throughflow Formulation 

Equations (2.11) through (2.14) can be further simplified by forming an explicit scheme that is not 
differential in nature. The computer code MTFLOW implements such a formulation [3]. 

The meridional flow speed, v m is obtained from a local streamtube conservation of mass. 



Q 

pA ( 2 7r 7* — BTe) 



(2.15) 



The streamwise momentum equation solved within MTFLOW is of the same form as Equation 

( 2 . 12 ). 



dP + pv m dv m 4 pvedve + Pd(AS) — pd(AW) — 0 (2.16) 

The change in work AW enters the formulation through the Euler turbine equation. 

A W = j Vd(rve) (2.17) 

If entropy losses or enthalpy additions due to heat release occur within the flow domain, they 
enter the formulation in Equation (2.16). Equation (2.16) is discretized in a finite volume sense as 
described in section 2.2.2 to preserve shock capturing. 

Instead of solving the differential energy equation, which would require large computational times 
for iterations, MTFLOW solves for the total enthalpy at any location in the flow. 

h 0 = h inl + ] -vf nl + AH + AW (2.18) 

The final result is a computationally efficient methodology for solving the inviscid equations of 
fluid motion, which allows for the effects of entropy, enthalpy, and rothalpy to enter the problem. 
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Inviscid And Solution Streamline Positions 




Dashed - Final Displaced Streamline Positions 
Solid - Elliptically-Generated Inviscid Streamlines 



Figure 2-2: The elliptically-generated streamlines and final streamline positions for the displaced 
body solution for Huang Submarine Body 1. Notice the overlap area between the displacement 
thickness and boundary layer thickness. 

2.3 IBLT Boundary Layer Representation 

While the fluid velocity and pressure gradients along the body are known from the inviscid streamline 
curvature solution, it remains to solve for the flow within the boundary layer. Drela’s formulation of 
the problem uses an integral boundary layer (IBLT) solver to find the boundary integral quantities 
knowing the inviscid edge velocity and pressure gradients along the solid boundary. Once the 
boundary layer flow is known, the body boundary is displaced by the boundary layer displacement 
thickness, and the inviscid outer flow is resolved. In this manner the boundary layer solution and 
outer potential flow solution are coupled together. 

2.3.1 IBLT Diffculties 

The difficulty for the propeller designer with the IBLT is that the IBLT solution gives the integral 
quantities of the boundary layer, such as momentum thickness and displacement thickness. Because 
the propeller design code requires the total inflow velocity over the entire surface of the blade, and 
portions of the blade extend through the boundary layer, the velocity profile of the fluid through 
the boundary layer has to be reconstructed. 

The second difficulty encountered is that the boundary layer displacement thickness, &* , is less 
than the actual boundary layer thickness, S. This is an issue, because the MTFLOW flow solver 
solution shows high inviscid velocities in the region between rf* and 6 , while in reality, the slower 
boundary layer fluid velocity exists. 

Both of these problems and solutions are highlighted in figure 2-2. 
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Chapter 3 

Coupling PBD and MTFLOW 



3.1 Process Overview 

The coupling between the streamline curvature code and the propeller blade design code follows that 
used by Kerwin [14] when he coupled the propeller blade design code with a Reynold’s- Averaged 
Navier-Stokes (RANS) code. The basis of Kerwin's coupling was the use of distributed body forces 
within the RANS code to represent the presence of the propeller. In the streamline curvature code, 
the propeller swirl, or angular momentum, is used to represent the presence of the propeller in the 
flowfield. 

In the present coupled analysis or design, one starts by using the propeller blade design code to 
predict the loading of the propeller based on an assumed flow field. The blade swirl is then transferred 
to the streamline curvature code where a throughflow solution is computed. To complete the cycle, 
a total inflow velocity field is extracted from the throughflow domain and returned to the propeller 
code. The induced velocities predicted by the previous run of the propeller code are subtracted from 
the total inflow velocity to arrive at a new effective inflow velocity field. This coupling methodology 
is shown in Figure 3-1. 

The coupling is automated through the use of small FORTRAN computer codes. The two 
computer programs shown in Table 3.1 were written to accomplish the coupling. The actual sequence 
of running computer programs to run a coupled propeller problem is shown in Figure 3-4. 



3.2 Obstacles to Coupling 

There were several obstacles to overcome in bringing about the PBD-MTFLOW union. 

1 . Converting propeller blade force quantities into a measure of energy addition to the surrounding 
fluid consistent with the streamline curvature basis of MTFLOW. 

2. Reconstructing a boundary layer velocity profile based upon integral boundary layer quantities. 



PBD2MT 


Converts the PBD output circumferential mean induced velocity to swirl. Write 
out the MTFLO input ascii file containing the swirl, and thickness, if present in the 
PBD run . 


BL2BODY 


Converts the MTFLOW output boundary layer data into velocity data (i.e . , bound- 
ary layer reconstruction), and writes out velocity data field for use by VELCON. 



Table 3.1: These two codes are the coupling routines written specifically to couple PBD and MT- 
FLOW. While no changes were made to any PBD code, the MTFLOW code was altered to output 
TECPLOT-formatted data files to affect the data transfer. 
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Propeller Blade Solver 




Streamline Curvature to 
Vortex Lattice Conversion 

O Boundary Layer Reconstruction 
Velocity Field Extraction 



♦ 



Streamline Curvature with 
Imposed Swirl 





Vortex Lattice Lifting Surface 
to 

Streamline Curvature Conversion 



O Blade Placement in Flowfield 
O Circulation to Swirl in Flowfield 



Figure 3-1: The coupling procedure between the propeller blade design code and the throughflow 
solver. 
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3. Extracting flow field velocities in a manner such that PBD internal data representation did 
not fail. 

4. Allowing for the severe radial redistribution of streamlines in the case of an open propeller. 

The solutions presented in the following chapters represent my proposed solution to the major 
obstacles. The final chapters, which present validation and results, show that perhaps there is some 
merit to these solutions. 



3.3 Non-Dimensionalization Issues 

The important flow quantities passed between the flow solver and PBD are the flow velocities, the 
propeller rotation rate, Q, and the propeller swirl, r&Vg. The flow velocities are used in the blade 
design code to enforce the kinematic boundary condition on the blade surface. The propeller rotation 
rate and swirl are used within the throughflow solver to represent the energy imparted to the fluid 
by the action of the propeller. 

3.3.1 Flow Velocities 

The magnitude of the velocity is always non-dimensionalized by the inflow speed. This is consistent 
between PBD and MTFLOW. Therefore, the velocities do not need to be altered when coupling the 
two codes, even in the case of contract ing or expanding, internal, or ducted flows, or in the case of 
internal flows where the domain inlet velocity is different from 1.0. 

3.3.2 Propeller Swirl 

The swirl is used in the Euler turbine equation to calculate the work done on the fluid by the rotor. 
Because of the length scale inherent in the term, the swirl must be scaled by the proper reference 
length scale when transposing swirl from PBD into MTFLOW. 

3.3.3 Propeller Rotation Rate 

Internal to PBD, the non-dimensional advance coefficient, J, is always given as follows. 

j = -h- 

nD 

1.0 2tt 

“2 R (2 

J=- 

fi 

V 8 is used to non-dimensionalize the velocities. V 9 is usually defined as the speed of advance 
of the vehicle. In the computation domain where the body is fixed and the flow field has non-zero 
velocity at infinity, V s is the value of the upstream velocity at infinity. PBD always assumes that 
V $ — 1.0. It is not correct to use V a to non-dimensionalize velocities to arrive at the advance 
coefficient. Notice that the value of the effective inflow velocity in the plane of the propeller swept 
volume, commonly called the apparent velocity V a , is not constant throughout the swept volume of 
the propeller. Hence, the use of Va would result in ambiguities between experimental and numerical 
predictions. 

It is common to non-dimensionalize empirical propeller tunnel test data with respect to the 
upstream velocity at infinity, V s . In internal flow cases, the upstream velocity at the inlet may not 
be equal to 1.0. Yet PBD assumes that V s is 1.0. Hence, if V s is not equal to 1.0, this fact must be 
reflected by altering the advance coefficient in PBD. This is the so-called effective advance coefficient . 
The use of an effective advance coefficient requires a change in the usually straightforward calculation 
of the rotation rate, fi. 



(3.1) 

(3.2) 

(3.3) 
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Je ff =^r (3.4) 

nD 

Q=i ^Vmle 1 (3.5) 

*^eff 

3.4 Open Propeller Issues 

Coupling PBD with MTFLOW to solve for an open propeller presents a challenge because the 
streamline positions evolve with the solution. To maintain radial equilibrium, the streamlines with 
high induced mass flow bunch near the centerline unless constrained by a duct. This draws in the 
lower velocity freestream streamlines to the tip region of the propeller. 

The input flow parameter distributions (of. say, swirl rAV # , and thickness) are specified with 
respect to either a fixed (c, r) grid location or a streamtube. Modelling a rotor (or impeller) presents 
the competing requirements that the swirl input by the rotor does not change physical location, 
while the shed vorticit.y (which I also refer to as convected swirl ) naturally stays within the same 
streamtube, which radially contracts and changes position as the streamtube solution evolves. If the 
shed vorticit.y does not remain aligned with the local flow, then the wake is not force free. 

The problem when using a streamtube solver is that in the case of an open propeller, the stream- 
tubes will constrict in the region of the propeller due to the increased axial velocity of the fluid. 
This contraction is shown in Figure 3-2. The solution is the use of an elongation factor as explained 
in Section 3.4.1. 



3.4.1 Solutions for Open Propellers 

The MTFLO program gives the user the opportunity to specify that input flow changes are fixed in 
z — v coordinates or in the streamline s — t coordinate system. Fixed in the s — t coordinate system 
means that the absolute position input flow quantity will evolve as the streamlines evolve with the 
solution. 

In the case of swirl input to the flow domain by a rotor, or removed by a stator, if is not obvious 
whether the input positions should be fixed or evolve. For while the position of the mechanical 
blades are fixed, the swirl convected downstream certainly is not. 

The solution is to let the swirl evolve with the flow. To counteract the contraction tendency, the 
radius of the output swirl should be linearly stretched from the hub upward. In other words, with a 
fixed hub radius, all blade radii are linearly scaled such that the tip radius is set to the elongation 
factor. Now, as MTFLOW contracts the streamlines towards the hub, the swirl from the blades 
is moved into the correct geometrical position! The elongation ratio, or stretching factor, is found 
by a trial and error process. The best way to watch this while running MTFLOW is to guess an 
elongation ratio, run MTFLOW, and then view the output swirl distribution. The swirl input from 
the tip of the propeller should be at the geometrically-scaled 1.0 radius value. The proper use of the 
elongation factor is shown in Figure 3-3. 



Since the radial position of the swirl is changed, the swirl input into MTFLOW must be changed, 
too, such that the strength of the circulation T is conserved. Because there is only one length scale 
involved, it is a linear rescaling. 

Notice from Figure 3-2 that not all streamline positions are shifted an equal amount radially. 
Therefore, a refinement upon the linear rescaling from hub to tip is to alter the radial position of 
each streamline based on the local streamline shift, vice the global shift. To actually find the shift 
in streamline position, this thesis implemented a procedure to examine the radial position of the 
circumferential mean tangential velocity along the trailing edge of the blade output by MTFLOW 
against the velocity distribution output by PBD. The difference in velocities at any radial positon 
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Tip Contraction 




Figure 3-2: Streamline contraction due to the presence of the propeller. Within MTFLOW if 
the input blade swirl is specified relative to streamline positions, vice physical (z,r) position, the 
blade swirl is unphvsically redistributed with the flow. The the blade design program (PBD) and 
MTFLOW are solving different flow problems. The surprising this is that the coupling still quickdv 
converges to a non-physical solution. 
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Figure 3-3: . 

Use of swirl elongation factor to counter streamtube contraction. Notice that the contour plot of 
swirl on the left hand picture is contracted such that the tip radius is less than 1.0. By increasing 
the elongation factor, the right hand side plot shows that the final tip swirl postion is nearly at the 
physical propeller radius of 1.0 



shows the radial redistribution of streamlines, and gives rise to a relaxation scheme to reconverge 
the streamline positions. 

The non-linear radial redistribution scheme is clearly superior for the case of a single blade row 
propulsor. However, it is not clear that this scheme will still work for multiple blade rows, where 
the swirl shed by the upstream blade rows must convect, and change the input quantities. While 
the relaxation scheme would certainly still work (the swirl input into the MTFLOW domain should 
be redistributed until it is not changed by the MTFLOW solution), it is not clear that convergence 
can be obtained as quickly, if at all. 

3.4.2 Swirl Convection 

The user is required to manually convect trailing swirl from the trailing edge of the blade row to the 
domain exit. While this is a well-defined problem for internal flow passages, where the streamline 
positions can not radically change, the open propeller ( i.e., external propulsor) presents unique 
challenges. From basic empirical observation as well as numerical simulation, it is known that the 
wake of an open propeller contracts. To maintain a force free wake, the swirl convected within the 
wake must contract with the streamlines. 

Several coupling methods were tried to overcome the difficulties with wake convection. 

1. Convect trailing wake with no radial contraction. This solution is the easiest to implement 
however it results in the physically unrealistic flow situation that the wake is no longer force 
free. This line of thought was not pursued to investigate the magnitude of the problem. 

2. Convect Wake Along previous iteration's streamlines. MTFLOW is used to solve the flow 
solution with the convected swirl fixed in (^,r). The convected swirl is placed downstream 
based on the previous iterations streamline positions. After the MTFLOW solution converges, 
the wake is re-aligned based on the updated streamlines, and the solution is re-converged. This 
wake re-alignment procedure is repeated two times before a new propeller blade solution is 
computed. 
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Run PBD 



Run PBD2MT 
Run MTFLO 

Run MTSOL 
Run BL2BODY 
Run VELCON 



Run PBD2MT 
Run MTFLO 

Run MTSOL 
Run BL2BODY 
Run VELCON 



Wake Alignment 
With Freestream 
Iterations 



Figure 3-4: The coupling procedure between the propeller blade design code and the throughflow 
solver for an open propellor case. The inner iteration loop is required to align the trailing vorticity 
with the freestream velocities. 



The final method works the best from a computational efficiency and accuracy standpoint. The 
iteration cycles are illustrated in Figure 3-4. 



In practice, after the initial MTFLOW solution has converged, the wake streamlines change very 
little, and it is practical to skip the inner wake alignment loop shown in Figure 3-4. 

3.4.3 Swirl Radial Redistribution 

Circulation shed off the propulsor blades is represented as swirl convecting along streamlines in 
MTFLOW. As the radius of the streamline about the machine's centerline changes, the strength 
of the swirl, rAV$, must change to maintain constant circulation, T. There is a direct inverse 
relationship. 



T oc / (rV$) 

This is an important consideration when convecting swirl through the meridional passage of a mixed 
flow machine and along the wake of a free propeller. If the swirl strength in the MTFLOW input 
were not altered to maintain constant circulation, the wake in MTFLOW would appear to have 
windmills removing energy in it as the radial position was altered. Recall that there can be no 
radial streamline redistribution due to wake vorticity, because there is no poewr being put into the 
fluid, as evidenced by the fact that Q equals zero in the trailing wake system. 
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./ 




Copy /Rotor/tflow.xxx ./rotor. tflow 
Copy ./Stator/tflow.xxx ./stator.tflow 
Run BUILDTFLOW 
Copy ./tllow.out ./tflow. xxx 
Run MTFLO 
Run MTSOL 

Copy VEUOIN.tev ./Rotor/VELJOlN.tec 
Copy VEUOIN.tec- ./Stator/VELJOIN.tec 


/ 


\ 


./Rotor/ 


./Stator/ 


Run VELCON 
Run PBD 14 
Run PBD2mt 


Run VELCON 
Run PBD14 
Run PBD2mt 







Figure 3-5: Multiple Blade Row Coupling Flowchart. This process is easily automated in a script 
file. 



3.5 Running the Coupling 

PBD2MT automatically assumes that the blade row described in the PBD output files is the only 
blade row that exists in the throughflow domain described in the MTFLOW walls. XXX file. If more 
than one blade row is present, PBD2MT must be run on each blade row, and then the BUILDT- 
FLOW program must be executed. This process is shown in figure 3-5. 



3.5.1 Coupling Admin FilerMTCOUPLE.INP 

The mtcouple.inp file is the only file which the user must create from scratch. It serves as an 
administration file. 



LINE 1: 3 

LINE 2: 0.00955 

LINE 3: 0.4211 

LINE 4: 9 

LINE 5: 0.2 

LINE 6: 1.0 

LINE 7: tf low. rotor 

LINE 8: walls. rotor 

LINE 9: rotor. pbd 

LINE 10: 1.000 



Reynolds number 

inlet mach number 

Vship used in PBD to calculate Js 

aft most point of body 

x location of LE tip 

r location of LE tip 

tflow file name 

walls file name 

pbdl4.4 input file name 

Elongation factor 



The input lines have the following uses within the code: 
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LINE 1 
LINE 2 
LINE 3 
LINE 4 

LINE 5 

LINE 6 

LINE 7 

LINE 8 
LINE 9 

LINE 10 



REYNOLD’S NUMBER. Reynold's number for the throughflow domain. If less 
than 10, BL2BODY code assumes an inviscid flow domain and no boundary layer 
reconstruction takes place. 

INLET MACH NUMBER. For incompressible flows, the mach number is defined 
as the velocity divided by the speed of sound. For water, M = x • See section 
A. 5 for further details. 

VSHIP. This is the upstream velocity at infinity. For coupling purposes it is the 
fluid velocity at the domain inlet. Remember that the advance coefficient, J , within 
PBI) is calculated using Vship. 

AFT MOST BODY POINT. Only used in the case of an external flow domain where 
the stern of the body terminates in the throughflow domain. Downstream of this 
point, BL2BODY will use its own wake model to reconstruct the boundary layer 
velocity profiles vice Swafford relations. 

X-LOCATION OF LE TIP. Given in throughflow coordinates, the axial location of 
the propeller leading edge tip - similar meaning as in VELCON. 

R-LOCATION OF LE TIP. Given in throughflow coordinates, the radial location 
of the propeller leading edge tip - similary meaning as in VELCON. 

TFLOW FILE NAME. Name of the MTFLOW tflow data file which will be output 
by PBD2MT. 

WALLS FILE NAME. Name of the MTFLOW walls geometry data file. 

PBD14.4 INPUT FILE NAME. Name of the PBD14 admin file. Used to read the 
advance coefficient and number of blades. 

Elongation factor as explained in Section 3.4.1. 
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Chapter 4 



PBD to MTFLOW Conversions 



4.1 Program Overview 

The coupling between PBD and MTFLOW occurs using the PBD2MT program. The file passing 
and interaction between PBD and PBD2MT is shown in Figure 4-1. The goal of the PBD2MT 
program is to produce an MTFLOW input ascii file, referred to as a tflow file in MTFLOW parlance, 
that contains a description of the swirl, entropy, thickness, and enthalpy input into the flow field 
by the propeller blades. Currently, PBD2MT handles swirl and thickness only. Entropy generation 
from the viscous drag on the blades could be added if desired, and enthalpy is not applicable for a 
propeller if the propeller forces are transferred via swirl. 

PBD2MT converts the PBD14-calculated circumferential mean propeller induced velocities to 
swirl, rV$ . The program also automatically checks for blade thickness and writes the blade thickness 
to the output file, as well. The user is only required to make a single administration input file to 
run the code. It would be a relatively simple matter to input the entropy generated by the blades 
based on the blade sectional drag into the meridional flow using the current framework. 



4.2 PBD2MT Input/Output Files 

The only specific input file the user must create is the coupling administrative file, called mtcou- 
ple.inp. The rest of the files are either created by the successful completion of a Mode 6 PBD 14.4 
run, or files required by MTFLOWL All of these files should be placed in the same directory where 
the PBD 14 output is written to ensure data is not lost or inadvertently overwritten. 



mtcouple.inp 


administrative file for PBD2MT program 
Lists flow paramters and input file names 


XXX.pbd 


PBD14.4 style input file 

Lists the Advance Coefficient, Js 


PBDOUT.CMF 


output file from blfrc.f with grid and control point 
(CP) positions as well as azimuthal blade thick- 
ness at CPs 


PBDOUT.CMV 


circumferential mean induced velocities at CPs 


WALLS.xxx 


MTSET style walls file 


restart.vel 


flowfield velocity description in tecplot format 
used for wake trajectory alignment 



Table 4.1: PBD2MT required input files 
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Figure 4-1: The interaction and file passing between PBD and PBD2MT. Notice that four external 
files are required to run the coupling: the PBD input file, the MTFLOW WALLS. XXX file, the 
fiowfield file ( usually called restart. vel), and the PBD2MT admin file, mtcouple.inp. The end 
product is the ascii MTFLOW input file, TFLOW.XXX 



tflow.XXX 


TFLOW file containing swirl and thickness distri- 
bution in the flow domain. Inputs directly into 
mtflo. 


p2mout.sgr 


The following output file is commented out in the 
code but could be used for consistency checks. 
Calculation for blade circulation at the trailing 
edge in same format as PBDOUT.SGR. 



Table 4.2: PBD2MT Output Files 



4.3 Program Flow 

Three main tasks are required to convert PBD- 14 output data to a form suitable for MTFLOW use. 

Task 1 read the walls. XXX file with the body geometry 
Task 2 Rescale the propeller geometry to the throughflow 
geometry scale 

Task 3 write a tflow.XXX file with the prescribed swirl, 
thickness, and entropy from the PBD output 

4.3.1 Task 1 : Read walls. xxx 

The walls. xxx file contains the data points of the solid boundaries. This same file is used by MTSET. 
PBD2MT reads the inlet and outlet x locations, and the solid boundary points. Since allocatable 
arrays are used in PBD2MT, the maximum number of points on any element is set by the variable 
MAXS1ZE in PBD2MT. If it is exceeded an error message is produced. 

For increased accuracy, the user should place many points along the hub and duct in the area 
of the blade, because PBD2MT uses linear interpolation to avoid spline interpolation overshoots. 
Specifically, a boundary point should be placed immediately upstream and downstream of the blade 
hub and tip to ensure a good interpolation. Also, the user should well-define any regions of large 
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curvature gradients in the hub and tip region of the domain. 

4.3.2 Task 2 : Geometry and Swirl 

The TFLOW.XXX file involves specifying the swirl r(AV^), thickness, and entropy distribution 
throughout the fluid domain. PBD2MT computes these quantites at the PBD blade control points, 
as explained below. The heart of this process is the aoutput.f subroutine. 



Blade Row Rotation Rate 

PBD2MT reads the PBD Administration file to get the advance coefficient and the number of blades. 
The blade row rotation rate is computed from Equation (4.1). 

a = P HD 

Jeff 

As explained in Section 3.3.2 Vs may not equal 1.0. Hence, the value for Vs must be included in 
the uncouple admin file. If J e jj > 10.0 then PBD2MT assumes a stator. 

Geometry 

PBD2MT reads the PBDOUT.CMV and PBDOUT.CMF files to get the following data: 

1. X,Y,Z location of lattice grid 

2. X,Y,Z location of control points 

3. (from .CMF) Azimuthal blade thickness at each control point 

4. (from .CMV) u,v,w circumferential mean induced velocity components at the control points 



Circulation to Swirl Conversion Equations 



The governing equation for the conversion of circulation to tangential swirl velocity is repeated here 
for clarity. 



f B (r) 



—2irrVg 

Z 



(4.2) 



Reversing the terms in Equation (4.2) gives the relationship for the tangential swirl in terms of 
circulation. 



-v; 



F B (r)Z 
2tt 



Equation (4.2) and Equation (4.3) are both required in the conversion process. 



(4.3) 



Circulation to Swirl Implementation 

There are five basic steps to convert circulation to swirl. 

1. PBD outputs the circumferential mean induced velocity u,v,w components at every control 
point on the blade in the PBDOUT.CMV file. A transformation to cylindrical coordinates 
yields the circumferential mean induced tangential velocity component. Now, the circulation 
at every control point is known through equation (4.2). 

2. The circulation T is scaled by the reference length common between PBD and the through- 
flow domain. This reference length is always the propeller radius in thr.oughflow coordinates, 
because in PBD this length is always 1.0. Circulation T does not usually need to be scaled by 
a reference velocity, since in both PBD and MTFLOW, all velocities are non-dimensionalized 
by the inlet velocity. 
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3. In the case of an open propeller the swirl may have to be elongated, as explained in Section 
3.4.1. To maintain constant T, if the radius is increased, the swirl rVs must be decreased. 

4. After non-dimensionalizing all arrays containing a length scale, the swirl, array name RVT m 
the source code, is then re-calculated using Equation (4.3). Finally, because MTFLOW is a 
compressible flow code, the isentropic compressibility factor must be applied to the swirl. 

r(AVe) = RVT(l.0 + 0.2A/ 2 ) (4.4) 

4.3.3 Task 3: Output to the tflow file 

The tflow. XXX file is broken down into blocks of data, which correspond roughly to span lines 
of constant percent chord length along the blade. In PBD14-speak, these are lines of constant n. 
Whenever there is an abrupt change in one of the specified flow variables, such as at the leading 
edge or trailing edge of the blade, a spline break is introduced by doubling that span line. 

One implementation key is that MTFLOW can only accept a set number of t values. The radii of 
the control points in PBD are determined by the local streamline vice a constant distance off of the 
centerline of the machine. The distributions must be output along constant spaced t coordinates. 
This requirement requires a cubic interpolation of the swirl from the PBD output radial position to 
the constant i coordinate. 



4.4 Coversion to MTFLOW (s.t) Coordinate System 

One key to the coupling process is the rescaling of the propeller swirl such that in MTFLOW s — t 
coordinates, the radial t coordinate goes from 0.0 at the hub to 1.0 at the edge of the domain. This 
means that if the user writes out a r coordinate in the tflow. XXX file, r, too, must go from 0.0 to 
its ultimate value. In general, the relationship between t and r is shown in Equation (4.5). 



U = 






o 



(4.5) 



The local streamline radius, r,-, is calculated by first dividing the flow channel between the hub 
and tip into equally spaced subdivisions, and then finding the non-dimensional distance of each 
subdivision from the hub. 



_ Rtip Rhub 

NJ 



where N I is number of radial positions 



iAr 

Rtip ~ Rhub 

where i is local radii index 



4.5 Important Usage Tips 

4.5.1 walls. xxx File 

If an axisymmetric element extends to the upstream and downstream extent of the fluid domain 
,such as in the case of the inlet and outlet of a waterjet, extend the geometry definition of the hub 
and casing slightly ahead of the XINL specified in the walls file. This will keep PBD2MT from 
having problems of finding the inlet radius. If this is not done, PBD2MT will write out a very 
strange tflow file. 
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4.6 Notes on the code PBD2MT 



4.6.1 Main program: pbd2mtv2.f 

1. Head the mtcouple.inp file. 

2. Read the WALLSNAME file. Converts the counterclockwise ordered data points into inlet-to- 
outlet. ordered boundary points. 

3. Read the pbd admin file to get number of blades and rotation rate. 

4. Read PBDOUT.CMF to get blade locations, CP location, azimuthal thickness. 

5. Read PBDOUT.CMV to get the TANGENTIAL induced swirl velocity. 

6. Read restart. vel to get streamline information. 

7. Rescale streamlines by SFAC. 

8. Convert length scales from PBD to MTFLOW scale. 

9. Call subroutine aoutput to write out tflow file. 

4.6.2 Subroutine aoutput.f 

This subroutine bears the responsibility for writing out the ascii tflow. XXX file for use in MTFLO. 
It has several components. 

1. Calculate the constant 1 and r values where data will be output. 

2. Write out the header lines of the tflow file. 

3. Write out a bunch of zeroes at the inlet. 

4. Write out zeroes at LE of blade. 

5. Double LE of blade and write out correct rotation rate. 

6. Loop down every blade span line and write out by 

(a) Spline x,r v$, and thickness with respect to r. 

(b) At every output r radius, interpolate the splines and write out values of x,r v$, and 
thickness. 

7. Double the specified value at the blade TE. 

8. Write out user specified number of convected wake stations between blade TE and domain 
outlet. 
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Chapter 5 



MTFLOW to PBD Conversions 



The main challenges in converting MTFLOW output back to PBD input was the reconstruction of 
the boundary layer velocity profiles, and matching these boundary layer profiles with the displaced 
body inviscid velocity flow field found by MTFLOW. The flow of input and output files is shown in 
Figure 5-1. 



5.1 BL2BODY Guide 

The following two sections present how to run BL2BODY and then theory behind the actions 
performed by BL2BODY. 

5.1.1 Running BL2BODY 

The input files required by BL2BODY are shown in Table 5.1. The output files created by BL2BODY 
are shown in Table 5.2. All files are in TECPLOT format. 



5.2 Velocities from MTFLOW to PBD 

PBD expects to receive the nominal wake from the flow solver. Internal logic within PBD extracts 
the effective wake. This requires an accurate representation of the velocity field in the region of the 
propeller, and downstream of the propeller. Therefore, BL2BODY must reconstruct a boundary 
layer velocity field through the IBLT displaced body region and attach it with the inviscid solution. 



5.3 Boundary Layer Reconstruction 

The Integral Boundary Laver solver within MTFLOW solves for the gross boundary layer quantities 
such as shape factor //, momentum thickness 0 , displacement thickness, 8 * , edge velocity, and shear 



1. mt couple. inp 
1. OUTBL.tec 



2. O RIG GRID, tec 

3. OUTVEL.tec 



Coupling administrative file. Fully described in section 3.5.1. 

Boundary Layer Quantities from MTFLOW: XSSI, H, DHSDHK, THETA, 
DSTR, UEDG, CTAU, TAUW, DISPq. Output by MTSOL. 

Original elliptic grid set by MTSET. Output by MTSET. 

Outer, inviscid, velocities from MTSOL solution. Output by MTSOL. 



Table 5.1: BL2BODY required input files. 
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VELCON7 style input file 



Figure 5-1: The interaction and file passing betweenMTFLOW and BL2BODY. The end product 
is the ascii, TECPLOT format file VELJOIN.tec with the velocity field. This file is formatted as a 
VELCON7 input file. 



OUTPROFILE.TEC 

BLZONE.TEC 

VELJOIN.TEC 


Boundary layer velocity profile as separate zones 
Boundary layer velocity profiles as one large zone 
Makes for easier viewing of contours, etc . . . 

Final flow field velocity with boundary layer and and inviscid flow joined 
together. Input file for VELCON71. 


dstarcompare.tec 

thetacompare.tec 


Debugging Analysis files 

MTFLOW displacment thickness 6 * along the body and BL2BODY final 
output 6 * 

MTFLOW momentum thickness 0 along the body and BL2BODY final 
output 0 



Table 5.2: BL2BODY Output Files. 
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stress coefficient. Because often times the root of the propeller blade extends through the boundary 
layer, it is necessary to recreate the fluid velocity profiles through the boundary layer. Indeed, in 
the full stern submarine concept, the entire propulsor is within the boundary layer. At the heart of 
the boundary layer relations is the closure relation. The current implementation uses the Swafford 
correlations as closure relations. 

5.3.1 Swafford Boundary Layer Profile 

Swafford presents an empirical correlation method relating the fluid velocity ratio across the 
boundary layer to the profile shape factor, H , and momentum thickness Reynold’s number based 
on edge velocity Re$ [23]. The methodology for genarating boundary layer velocity profiles is well 
documented in the previously cited reference. These velocity profiles are valid for attached and 
separated turbulent flows. 

The velocity profile within the boundary layer is first computed in local body normal £- rj coordi- 
nates using Swafford’s correlation. Values for the local shape factor, //, and Ree are provided from 
the MTFLOW solution. Swafford’s procedure yields the total velocity magnitude. This velocity 
profile is checked to ensure that it matches the displacement thicknes, rf*, output by the MTFLOW 
solution. If it does not, the profile shape is linearly adjusted until the displacement thickness matches 
that of the MTFLOW solution. 

5.3.2 Boundary Layer Wake Model 

Analytically resolving the boundary layer velocity profile in the wake near the trailing edge of a 
body has been a great source of difficulty for those researchers not pursuing a volume representation 
of the fluid domain. A similarity solution to the wake profile is valid in the laminar case only after 
three body, or chord, lengths downstream of the body. This fact leads most researchers to deduce 
an interpolation scheme from the body boundary layer, through the near wake area, to the far wake 
(3 body lengths downstream) where analytically-derived, empirical relations exist. 

For this research, the analytical approach derived by White was followed [25]. White shows that 
the velocity defect within a laminar axisym metric wake is given by: 



Ul _ U 0 LL U 0 r\ 
Uo ° ar eX ^ 4xi> ' 



(5.1) 



If equation (5.1) is applied to a turbulent wake, velocity ratios greater than 1.0 result. To avoid 
this situation, this research modified equation (5.1) such that a Gaussian distribution of velocity 
across the wake was maintained: 
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(5.2) 



where Xte 


X location of body TE 


/ ar 


X location where far wake approx, valid 


r 


local radius 


Q 


adjustment factor 



The adjustment factor a is adjusted such that the momentum thickness of the resulting profile 
shape is equal to the momentum thickness given by the MTFLOW soltuion. 
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5.3.3 Fitting Boundary Layer Velocity Field to Flow Field 

Using either the straight Swafford relations or the modified wake relations, a velocity is profile is 
found through the boundary layer. It now remains to body fit the £- 7 / coordinates to the body 
boundary. Since the boundary layer velocity profile is found in a normal coordinate system, it would 
seem natural to merely rotate the £- 7 / coordinate system such that the 77 normal direction vector 
matched the vector normal to the body. This method is incorrect, however, because the boundary 
layer velocity vectors do not conform to the flow streamlines. This method also suffers from the 
shortcoming that along concave body surfaces, the boundary layer profiles will intersect at a distance 
equal to the local radius of curvature off the body. 

A second method is to keep the velocity vectors normal to the body, but lay the velocity vectors 
along a gridline semi-normal to the body. This method prevents boundary layer velocity profiles 
from crossing, but does not ensure that the boundary layer velocity vectors follow streamlines. 

A good compromise, then, is to lay the velocity vectors along an elliptically-generated gridline 
semi-normal to the body, and force the velocity vectors of the boundary layer flow to exactly conform 
to the inviscid streamlines. Near the body, where the inviscid streamlines are displaced off the body 
by the displacement thickness, the streamlines from the original elliptically-generated grid are used 
the reader will recall that the grid points in an elliptically generated grid numerically satisfy 
Laplace’s equation, and, thus, are inviscid streamlines. This compromise method gives very good 
results. 

5.3.4 Joining Boundary Layer Velocity Vectors and Outer, Inviscid Ve- 
locity Vectors 

MTFLOW solves for the inviscid fluid velocity outside of the displacement thickness, S * . S* is always 
less than the boundary layer thickness. There is an overlap region , then, between the boundary layer 
velocity profile, which extends from the body to a distance 8 off the body, and the solved, inviscid 
velocity flowfield, which extends from a distance off the body to the edge of the fluid domain. 
The velocity within the overlap region was taken as the boundary layer velocity for 90 percent of the 
boundary layer region - arbitrary, but reasonable. To correct for any interpolation problems, the 
solved-for boundary layer velocity is adjusted by a constant factor such that the velocity at the edge 
of the boundary layer matches the inviscid velocity within a streamline near that point. Since this 
linear scaling alters the displacement thickness and momentum thickness of the boundary layer, the 
boundary layers must be adjusted to match the outer inviscid velocity and the momentum thickness. 
The final 6* of the velocity profile closely matches the MTFLOW-derived 8 * . 

Swirl in the Boundary Layer 

Because the tangential swirl velocity is not carried in the boundary layer quantities, the inviscid 
swirl must be smeared across the boundary layer to give good results. If the swirl is not represented 
correctly, the blade section within the boundary layer will have the wrong pitch angle. 

The tangential velocity is interpolated across the boundary layer. The tangential velocity along 
the boundary is considered 50 % of the tangential velocity at the inviscid edge of the boundary layer. 
I 11 between a linear interploation is done to transition for the value at the wall, to the value at the 
inviscid edge 



5.4 Notes on the code BL2BODY 

The code is written in FORTRAN 90. Where possible, allocat.able arrays are used and destroyed 
within individual subroutines. IMPLICIT NONE is used to ease debugging. All data is passed 
between subroutines as arguments. There are no common arrays. Where parameters are hard-wired 
into the code, the offending subroutine header is annotated, as well as comments placed in areas of 
the code affected. 
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5.4.1 Main Program : bl2body.f 

The file bl2body.f controls the sequence of events required to convert MTFLOW velocity and bound- 
ary layer data into a flowfield suitable for use by PBD-14. 

1. Read the administration file, mtcouple.inp 

2. If the Reynold’s number is less than 10.0 BL2BODY assumes that this is an inviscid flow case 
and calls outvel2veljoin. 

3. Read the boundary layer data file called OUTBL.tec which is output by the modified MT- 
FLOW lo.f subroutine. 

5.4.2 Subroutine sprof.f 

This subroutine computes the velocity at any point within the boundary layer based on the mo- 
mentum thickness-based Reynold’s number Reg and boundary layer shape factor H . The routine is 
attributable to Swafford [23]. If the routine has any problems matching the correlations, usually due 
to trying to evaluate the velocity in a region of laminar flow, the subroutine automatically switches 
over to a I th Power Law velocity profile calculation. 

5.4.3 Subroutine joiner.f 

The boundary layer velocity profile data previously found is passed in. Because the boundary layer 
data is already in the correct x,y position it now must be joined with the inviscid velocity solution 
from MTFLOW. The inviscid solution is written out by the modified MTFLOW subroutine lo.f in 
the file OUTVEL.tec. 

The steps accomplished in joiner.f : 

1. Read in the OUTVEL.tec invsicid velocity data. 

2. If there is fluid along the centerline ahead of the axisymmetric, centerline body (e.g.. in the 
case of a shaft with a propeller on it which extends from the downstream portion of the flow 
domain into the flow domain) interpolate the inviscid velocities and write to the output file. 

3. Shift the boundary layer velocity data such that the velocity is given at the cell centers, vice 
nodes 

4. For each vertical line of boundary layer data 

(a) Find the inviscid velocity point closest to the end of this line 

(b) Consider this point the inviscid edge velocity 

(c) Linearly alter the boundary layer profile previously found to match this inviscid edge 
velocity 

5. Calculate the new boundary layer displacement and momentum thickness for the user to 
visually compare with the MTFLOW solution for the displacement and momentum thicknesses. 

5.4.4 Subroutine velconout.f 

The only purpose of this subroutine is to take the joined velocity file written out by joiner.f which 
is j-ordered data, and turn it into z-ordered data. This subroutine: 

1. Open and read VELJOIN.tec 

2. Over- write VELJOIN.tec in i- ordered data format 
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Chapter 6 

Validation 



Validation of the MTFLOW-PBD coupling was carried out in several stages, following the U build 
a little, test a little” philosophy. These tests built from the testing of individual components, and 
progressed through the validation against empiricial test cases. 

Because the two underlying computer codes are well- validated, it was not necassary to verify the 
internal accuracy of the codes. Rather, logical consistency between the codes was required. This 
was done by using simple propeller models, such as one with zero thickness, inviscid blades, and 
showing consistency between inviscid theory and MTFLOW. 

The coupling technique was validated in several stages. Starting from validation of the basic 
flow solver abilities, each successive level of validation represents an added level of complexity. By 
the last validation, a fully viscous internal flow case is tested. Table 6.1 shows the validation tests 
conducted. 



6.1 Validation: Nominal Wake 

As previously mentioned, if there is no vorticity present in the inflow (z.e., an inviscid inflow with no 
radial velocity gradients), there is no vorticity to interact with the propeller, and no inflow vorticity 
to redistribute. Hence, the effective inflow must equal the nominal inflow. This fact is used in 
validating the coupling between the propeller design code and the throughflow solver to ensure that 
the same forces, or power in this case, are modelled within the throughflow solver and the propeller 
code. 

Modelling a propeller on an infinitely long hub of constant radius removes any hub and duct 
effects from the circumferential meridional flow. Thus, we can assume a fluid inflow velocity which 



Open Propeller, inviscid 


To validate the ability of the two codes to couple 
and converge the nominal and effective wakes. 


No propeller, viscid 


To validate the ability of MTFLOW to accurately 
predict low mach number, high Reynold's number 
flow situations. 


Open Propeller, viscid 


To validate the abilityof the two codes to couple 
and converge to experimental results 


Ducted propulsor, inviscid 


To validate the ability of the two codes to recreate 
empirical ducted propulsor results. 


Ducted propulsor, viscid 


To validate the boundary layer routines and cou- 
pling technique. 



Table 6.1: Validation of the coupling technique involved increasing levels of fidelity and complexity. 
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is constant across all radii. This is the nominal inflow. After the PBD-MTFLOW coupling has 
converged (usually after about 6 interations), we can extract the effective velocity from the total 
flowfield. If this effective velocity from the coupled solution matches the nominal velocity of 1.0, we 
will have shown that the codes couple to the correct answer. 

The propeller code solves for the blade circulation given the constant inflow. The circulation 
is input into the throughflow solver as swirl. The throughflow solver solves the equations of fluid 
motion. In so doing, the fluid velocity is increased due to the presence of the propeller, streamtube 
contraction occurs across all propeller radii, and wake contraction becomes evident downstream 
of the propeller. Once the flow solution is converged, the meridional velocity field in the region 
of the propeller is extracted, and converted for input into the propeller design code. Those are 
the steps which comprise one coupling iteration. Now, the propeller code re-solves for the blade 
circulation distribution with this new inflow. This process takes approximately six coupling iterations 
to produce converged results. Results are considered converged when the circulation distribution 
on the propeller does not change with successive throughflow coupling iterations. Conversely, the 
throughflow flowfield velocity solution does not change with successive couplings. 

6.1.1 KA 455 Rotor 

A Kaplan KA 455 rotor was used for this validation. The nearly radial projected blade shape reduced 
problems with MTFLOW splining. Figure 6-1 shows the results for a Kaplan Series 4-55 rotor with 
no duct placed in a nominal inflow of unit velocity. Notice the over the majority of the blade, the 
effective velocity is with 1 percent of the nominal velocity. The missed regions in the hub indicate 
that there may still be some interpolation problems which need to be smoothed out near the hub. 
This coupling was performed using the non-linear streamline radial redistribution methodology. 

6.1.2 Nominal And Effective Calculation on the Fly 

Since the nominal and effective inflow velocity fields must be equal for the propeller in an inviscid 
flow, there is another sure-fire method to test any new coupling routine. This is to compare the radial 
load distribution predicted by the coupled codes with the radial load distribution predicted with the 
propeller in the nominal inflow field, with the caveat that the wake trajectories from the current 
meridional flow solution are used to ensure that wake convection between the two calculations is 
exactly the same. If the solution is not converging rapidly this method very quickly eliminates the 
wake sheet geometries as the culprit, since the same wake sheets are used, now, in both calculation. 

One key is to get the same trailing wake sheet geometries between the effective calculation and 
the coupled calculation. This involves specifying the local wake convection axial and tangential 
induced velocity components in the PBD Mode 4 input file[ 1 6] . The radial contraction portion is 
accounted for by the fact that PBD evolves the wake along the input streamlines. To find the 
induced velocity component in the wake, the circumferential mean solution from the throughflow 
solver is used. Since the upstream flowfield is uniform with strength 1.0, any difference in velocity 
at any point from 1.0 must be due to an induced velocity component. Now that the circumferential 
induced velocity is known in the wake, the wake is thus aligned correctly. 

As a further refinement, the position of the streamline positions from the throughflow solution 
are also used. That is, the streamline positions from the VELCON output file are used. The velocity 
along these streamline is set to axial, uniform inflow. Now the user is assured that the vortex lattice 
geometry used in the coupled and uncouple, stand-alone effective wake calculation are the same. 

Once the radial circulation distribution is found from the uniform effective flowfield, the coupling 
is run. When the solution converges, since the nominal wake must equal the effective wake the 
propeller code should give the exact same blade loadings. The results for this test, again performed 
with a Kaplan KA 4-55 rotor, sails duct, are shown in Figure 6-2. The results from this analysis 
show a completely converged solution. 

In performing this validation, the required radial redistribution of swirl problem was highlighted. 
It became apparent that the linear redistribution of radii would not produce the correct results. 
Therefore, the radial posit ion of each constant span line of swirl input into MTFLOW was altered by 
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Effective V elocny o-n B lade after M TP LOW C oup lin g 
Invieeid, 2cro Ihiokncee 




Figure 6-1: . 

Given a nominal axial inflow velocity of 1.0, the axial effective inflow velocity shown in the contour 
plot should equal 1.0 along the entire blade. The small portions of the blade not equal to 1.0 , level 
3 on the chart, are due to inaccuracies in the present streamline adjustment methodology. 
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Nominal and Effective (Converged) Blade Loading 




Figure 6-2: . 

The PBD Mode 4 calculation was done with a 19x19 vortex lattice. This is the converged grid 
size. Since there is no vorticity ( i.e no boundary layer) redistribution of the nominal wake, the 
predicted radial loading distribution from the effective wake calculation (z.e., PBD Mode 4 run) and 
the coupled MTFLOW-PBD are identical. 



examining the difference between the circumferential mean tangential velocity profile at the trailing 
edge of the blade that was output from PBD and the circumferential mean tangential velocity profile 
at the trailing edge of the blade that was output from MTFLOW. The throughflow code can not 
alter the tangential circumferential mean velocity . To do so would require adding or subtracting 
power from the fluid. Therefore, the radial distribution of tangential velocity output by the propeller 
code must be the same as that within the throughflow solver. To ensure consistency, they should be 
identical. In practice, however, if they are identical, too much swirl is introduced into MTFLOW, 
and the blade overpowers the flow. Hence, the circumferential mean tangential velocity profile output 
by MTFLOW should be about 3-5 percent lower than the PBD input. A sample comparison velocity 
profile is shown in Figure 6-3. 



6.2 No Propeller, Viscous Validation 

The boundary layer routines themselves were validated using tests listed in Table 6.2. These are 
all essentially stand-alone tests and independent of the flow solver and the propeller code. As they 
represent more of a code debugging effort than a theoretical development, the resulting plots are 
not shown. 

6.2.1 Huang Body of Revolution Nominal Wake 

The Huang body of revolution is an axisymmetric, unappended subsea notional test vehicle tested 
extensively throughout the 1970’s and 1980’s [9], [8] [7], [6]. There were a total of three models built, 
which differed in the after body buttock lines. The primary use of these models was to investigate 
the boundary layer characteristics in the stern region of the model with the three different stern 
geometries, with and without a propeller running. The model geometry for the so-called Huang 
Body 1 is shown in Figure 6-5. 
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Vtheta at Blade Trailing Edge Output From PBD and MTFLOW 




Figure 6-3: The circumferential mean tangential velocity at the blade trailing edge position output 
from PBD and MTFLOW. They are not equal because the blade circulation extends beyond the 
blade tip, overpowering the fluid. Hence the overall strength of the swirl (circulation) was decreased 
to account for this effect 



1 


Reconstruction 
of Swafford 

Profiles 


The closure relations used in the MTFLOW for- 
mulation are provided by the Swafford profiles. 
The routines were required to exactly reproduce 
the Swafford profile for the test conditions indi- 
cated in [23]. 


2 


Comparison to 
k th Power Law 


Using the same inputs for the boundary layer, a 
profile was created using the i Power Law and 
the Swafford routines. 



Table 6.2: Boundary layer velocity profile validation tests. 
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Streamline Radial Repositioning 




Figure 6-4: This plot highlights the non-uniform repositioning of the streamlines required to perform 
the nominal wake validation test successfully. 
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Huang Body 1 

Length / Body Diameter 10.9735 
Length / Prop Diameter 20.1 35 



Model Characteristics 

Body Length 21.947 

Body Radius 1.000 

Propeller Radius 0.545 

Propeller Leading Edge Tip 21.630 



Boundary Layer Velocity Profile Extractions 




Figure 6-5: Huang Submarine Body with Stern 1. 



Measurements carried out on Huang Body 1 included measurement of pressure along various sta- 
tions, velocity profiles at multiple stations along the body with and without a propeller, and detailed 
measurements of boundary layer quantities necessary to refine boundary layer closure relations. 

6.2.2 MTFLOW as a Viscous Calculator 

While this thesis is concerned with the interaction of a propeller and free stream vorticity, MTFLOW 
was first validated as a suitable tool for analyzing the flow around an external body with no propulsor. 
As such, numeric tests were conducted on the unappended Huang Body 1 model to predict the 
nominal inflow velocity in the stern region of the body. The velocity profiles calculated from the 
present method with the IBLT and the Swafford reconstruction algorithm were compared with the 
results from two different RANS codes [24] and the empirical data, as tested by Huang [9]. The 
velocity cut is taken very near the stern of the body as shown in Figure 6-5. The comparison is 
shown in Figure 6-6. The discrepancies between the three numeric codes and the empirical data 
suggest that further refinements are possible in the realm of boundary layer modelling. 

6.2.3 V h Power Law Boundary Layer Comparison 

An inquiry was made into the general utility of using the complex IBLT formulation by comparing the 
results of the IBLT with a simple \ power law calculation of momentum thickness and displacment 

thickness along the Huang Body 1. Figure 6-7 shows that the I th power law calculation is very poor 
within the stern region, and any velocity predictions based upon this result would suffer accordingly. 
Also, a resistance prediction of the body based upon an analysis of the downstream momentum 
thickness would lead to an inaccurate thrust requirement for the propulsor. 



The difference in predictions on the aft part of the body show the influence of neglecting the 
pressure gradient due to curvature in the streamwise direction, and the crossflow pressure gradients 
arising due to transverse curvature. The resulting drag, computed from the downstream momentum 
thickness, is an order of magnitude different. 

6.2.4 Boundary Layer Velocity Extraction Validation 

To test the boundary layer velocity extraction routines, boundary layer velocity profiles were cal- 
culated using the Swafford routine and a simple j th power law prediction. In fact, White [25] 

reccommends use of the power law for general applications. In the ~ th power law there is a 
simple relationship between velocity and distance off the boundary. 



50 



Comparison of Body 1 Nominal Inflow at Propeller Plane 
RANS, MTFLOW, EXPERIMENT 




Figure 6-6: Comparison of MTFLOW, RANS, and Experiment Velocity Nominal Wake Profiles. 
Notice how well MTFLOW predicts the results in comparison to some of the RANS results, which 
take orders of magnitude longer to run. 
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Figure 6-7: Comparison of MTFLOW IBLT predicitons and I th power law prediction for 6* and 0 
along Huang Body 1 



There is also a simple relationship between the boundary layer thickness and the displacement 
thickness [18]. 



S = 7.987 6* (6.2) 

Figure 6-8 shows representative velocity profiles predicted by the Swafford profile relations and 
the i Power Law for the same input parameters of S* and u e d ge . The body used is the Huang Body 
1. The boundary layer quantities 8* and u € d ge were taken from the MTFLOW solution. Notice that 
far upstream of the propeller, where there are low pressure gradients, the A Power Law prediction 
is in good agreement with the Swafford correlation. Near the propeller plane at the aft end of the 
Huang Body 1 model, however, there are large discrepancies, especially within the model propeller 
radius of 0.545. 



6.3 Open Propeller, Viscous Validation 

The best way to test the performance of the present coupling methodology is to validate against 
experimental data. Two test cases were chosen to validate the coupling procedure. 

The first validation is the the Huang Body 1 tests with propeller 4577 operating at a constant 
J value. This represents an extremely hard case, because the propeller operates entirely within the 
boundary layer of the body. This requires a robust boundary layer method. 
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x/L = 0.65 







Figure 6-8: . 

Comparison of Swafford Profile and \ th Power Law Boundary Layer Velocity Profiles. The location 
of the stations are shown in Figure 6-5. The same boundary layer quantities predicted by the IBLT 
were input into a velocity reconstruction routine based upon the Swafford data, and one based upon 
the 4 Power Law. 



53 



J = 1 .25 


Experiment 


Current 

Technique 


Ax 

IOA'q 


0.2276 

unk 


0.2164 

0.4685 



Table 6.3: Huang Body 1 with Propeller 4557: experimental and numerical results. 



The second test case is the 1998 ITTC Propeller Workshop conference test case of the three 
bladed propeller 4119 on a downstream driven shaft. This case was chosen to allow comparison 
with the validation performed by ITTC. This canned test case really tests the interaction of the 
propeller and the inviscid flow, since the boundary layer is confined to a very thin region close to 
the shaft, and presents an opportunity to see how the coupling procedure performs under real world 
conditions. 

6.3.1 Huang Body 1 with Propeller 4577 

The stern section of Huang Body 1 with the projection of propeller 4577 is shown in figure 6-9. The 
gridlines in the Figure are gridlines, not flow streamlines, or finite volume cells. Figure 6-9 highlights 
the complexity of the propeller in the computational domain due to the non-cylindrical propeller 
sections, lion-radial propeller generator line, and non radially uniform inflow field. 

The particulars of propeller 4577 are shown below. 



Number of Blades 
Section Meanline 
Section Thickness 
Rake Angle (deg) 
Skew (deg) 



NACA a=0.8 

DTNSRDC modified NACA 66 

6.964 

30 



r 

Kp _ __ 


Chord 

Dp 


t 

c 


P 

D p 


1m. 

c 


0.2106 


0.171 


0.235 


0.823 


0.0014 


0.3000 


0.177 


0.209 


0.980 


0.0175 


0.4000 


0.182 


0.182 


1.151 


0.0288 


0.5000 


0.185 


0.158 


1.243 


0.0337 


0.6000 


0.185 


0.135 


1.264 


0.0341 


0.7000 


0.180 


0.116 


1.248 


0.0311 


0.8000 


0.164 


0.0995 


1.206 


0.0246 


0.9000 


0.132 


0.0875 


1.157 


0.0141 


1.0000 


0.069 


0.0813 


1.108 


0.0000 



The results from the coupling closely match experimental results, as shown in Table 6.3. 

6.4 Propeller 4119: ITTC 1998 Tests 

In Spring 1998, the International Towing Tank Committee (ITTC) held a propeller workshop for 
RANS codes and panel methods. The premise for the workshop was to validate various numerical 
codes against empirical test data for propeller DTMB 4119 tested in the Naval Surface Warfare 
Center, Carderock Division 24 inch water tunnel [10]. 

The test parameters for a single advance coefficient are shown below. 
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Figure 6-9: Huang Body 1 Stern Profile with Propeller 4577 



J = ^n 


0.83333 


V , advance velocity 


8-33^7 


??, rotation speed 


10 rps 


D 


12 in 


R 


6 in 


V »nie< 

r 


1.0093 



Note that the inflow speed at the test section inlet is not 1.0 . This required calculating an 
effective advance coefficient as described in Section 3.3.3. 

The following flow parameters were used within MTFLOW. 



Inlet Mach Number 


0.001488 


Reynold’s Number 


432000 


Far Field Type 


3 



The MTFLOW grid is shown in Figure 6-10. Grid is actually a misnomer, since the gridlines 
are actually the streamlines. As such, the initial grid represents the initial guess at the streamline 
positions. Streamlines were bunched near the shaft only to resolve the potential interference between 
the leading edge separation bubble and the propeller blade. The truncated grid domain was used to 
overcome MTFLOW spline interpolation problems. 
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MTFLOW Grid 

ITTC Downstream Driven Shaft with 4119 



Freestream Jet Boundary at R=2.00 




Figure 6-10: Propeller 4119 in ITTC Configuration 



J 


Exp 

A 7 1 


Present. 

A 7 


Exp 

IOA'q 


Present. 

IOA'q 


0.700 


0.200 


0.194 


0.360 


0.347 


0.833 


0.146 


0.141 


0.280 


0.269 


0.900 


0.120 


0.118 


0.239 


0.237 


1.100 


0.034 


0.045 


0.106 


0.131 



Table 6.4: Experimental and numerical propeller curve test results for 4119. 

6.4.1 Inviscid Comparison to Experimental Data 

An inviscid coupling was performed to validate the overall coupling methodology and investigate 
the magnitude of the radial streamline redistribution problem. 

The published experimental results for the Workshop test case and the numerical results from 
the inviscid PBD-MTFLOW coupling are shown in Table 6 . 4 . 

For this series of test, the linear redistribution of radii was used, vice the non-linear redistribution 
used in the KA455 test. These results show that the simplfied, linear redistribution is satisfactory 
to reproduce experimental results at lower J values. At high J the linear model results in too much 
loading at the tip, which results in too high AY and too high I\q. Of note, the calculation only 
took 2 hours to converge all data points. 



6.4.2 Problems and Solutions 

The biggest problem was that MTFLOW had a lot of problems creating a two dimensional spline 
representation of the input swirl. A sample spline interpolation problem is shown in Figure 6-12. 
The final solution was to alter the gridding until a “happy medium” was found between the PBD 
output swirl distribution and the MTFLOW grid. There did not appear to be a rational, consistent 
method to resolve the splining problems. 
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KT, 10KQ 



ITTC 1998 Propeller Workshop Propeller 41 19 




J 



Figure 6-11: comparison of MTFLOW-PBD with P4119 ITTC Test. 




Figure 6-12: This is a view within MTFLOW of a contour of the swirl input from PBD into 
MTFLOW. Notice the large jagged peak on the hub radii. This is the symptom of the spline 
interpolation problem. When the MTFLOW spline routine fails, it puts a 0 data point in the field. 
Thus, the flow field responds as if there were zero swirl in the flow at that point - a porous propeller. 
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Chapter 7 



MTFLOW-PBD 14 Fan Design 
Mode 



To test the robustness of the coupling methodology, the PBD-MTFLOW coupling procedure was 
applied to the analysis, re-design, and subsequent analysis, of a ducted fan unit. The design of a 
ducted fan presents several technical challenges: 

1. The high contraction ratio from the inlet to the fan section and vertical walls present in 
the transition section present difficulties to the flow solver. Actually, within PBD 14.4 it is 
impossible to extend the hub and duct vortex lattice system across a vertical boundary since 
all velocity data is interpolated via a cubic spline with respect to the axial location. On a 
vertical boundary, there are are multiple flow values at the same axial location. This situation 
is termed a multi-value solution and can not be resolved via the interpolation methodology as 
currently implemented within PBD 14.4. 

2. The very small radius of curvature duct lip just upstream of the leading edge makes the 
enforcement of the kinematic boundary condition within PBD on the duct difficult. 

3. Low inlet speeds can result in negative effective inflow velocities. Thus, the fan blades are 
fighting to overcome their own negative effective inflow. 



The particular fan unit design was motivated by a design competition sponsored by a local fan 
design firm [21]. The sponsor set the design rotation speed, flow rate, hub ratio, and duct geometry. 
All length scales were scaled such that the leading edge tip radius of the fan was 1.0 to be consistent 
with PBD scaling. 



7.1 Fan Radiator System Geometry 

The fan was designed to operate as part of a radiator unit, with the radiator upstream of the fan. 
The geometry is shown in Figure 7-1. This geometry was modelled in MTFLOW with the exception 



Rotation Speed 


2300 rpm 


Mass Flow Rate 


0.419 


Hub Radius 


0.40 


Tip Radius 


1.00 


Length Scale 


1.46 



Table 7.1: Fan Design Operating Parameters. 
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radius curvature = 0.015 




Figure 7-1: The geometry of the radiator fan system. This geometry was modelled in MTFLOW 
with the exception that the downstream hub and duct lip was modelled as a straight horizontal wall. 
Also, the radiator was not modelled. 



that the downstream hub and duct lip was modelled as a straight horizontal wall extending to 
infinity. This was done only to conserve computational expense. Physically, the high flow region 
around those lips is located too far downstream of the fan to interact with the fan blades. Also 
the radiator was not modelled. The effect of the radiator is to dump enthalpy into the fluid, and 
straighten the flow. So, the assumption is that the radiator is at infinity, and produces no change 
in enthalpy. 

The engine designer desires the shortest length fan possible to reduce engine “stack" length. 
In fan terms, we want the lowest possible pitch angles, such that fan stack length is not reduced 
by reducing the chord lengths of the fan blades. A reduction in chord lengths is undesirable for 
improved off-design performance. 

The final MTFLOW grid used is shown in Figure 7-2. The long downstream run is required 
by PBD to grow the trailing wake system. Final solution results showed little change in the flow 
characteristics downstream of one propeller radius. The effects of modelling a shorter wake within 
PBD were not investigated. The upstream length of the domain is set by the requirement that the 
inviscid streamlines correctly match the constant inflow assumption. 



7.2 Initial Blade Design 

The initial blade design was carried out by Dr. Robert VanHouten using his own two-dimensional 
blade section design code. Due to the internal geometry limitations of the code used, the original 
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MTFLOW Domain with Meridional Blade Superimposed 




Figure 7-2: The flow grid used within MTFLOW. The long downstream extent is required by PBD 
to grow the trailing wake system. The upstream length of the domain is set by the requirement that 
the inviscid streamlines correctly match the constant inflow assumption. 



fan design consists of seven radial blades. The initial fan design by Dr. Van Houten is shown in 
Figure 7-3. 



7.3 Radial Blade Analysis 

The initial radial blade was analyzed using PBD- 14.4 with a vortex-lattice hub and duct repre- 
sentation, coupled with MTFLOW. The blades were discretized using a vortex lattice mesh of 225 
elements. Similarly, the hub and duct were discretized upstream and downstream of the propeller. 
The hub and duct vortex lattice systems were not extended over the entire vertical sections upstream 
of the fan due to internal numerical limitations noted in section 7. 

The result of this inviscid analysis was a solved-for radial circulation distribution along the 
span of the blade. A comparison between the radial circulation distribution predicted by PBD14.4- 
MTFLOW and the original design circulation distribution along the span of the blades is shown in 
Figure 7-4. 



The final converged flow solution streamline mesh and axial velocity contours are shown in Figure 
7-5. The extremely high inviscid meridional velocity at the duct lip near the tip of the blade can lead 
to instabilities in converging the flow solution. This high velocity highlights one of the limitations 
of the current coupling technique. While the throughflow is for the circumeferential mean flow, 
the blade experiences a different local blade velocity, due to the three dimensional nature of the 
blade vortex sheets. While it is well-known that the three dimensional unsteady harmonics of the 
propeller solution attenuate rapidly away from the propeller, in the flow problem considered here, 
it is not obvious that there would be no interaction between the local, three dimensional, unsteady 
harmonic solution and the meridional flow which, due to the high velocity gradients in this region, 
may produce undesired loadings or blade fluctuations. 



7.4 40° Back Skew Blade Design 

The first design iteration was to design a fan blade with the same circulation distribution as the 
radial fan blade, but with the addition of 40° of backskew to the blade tip. The back skew serves 
two purposes: 
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2 




Figure 7-3: The original radial fan with cut-away duct. The trailing wake vortex system is shown 
emanating only from one blade to maintain clarity. 

Reduce non-circumferential mean effects . Quite often the propeller blade operates in a flow 
domain that is non-uniform. Usually this is due to an upstream obstruction or contorted flow 
passage. As the blade passes through this region of fluid with a different velocity, the apparent 
angle of attack of the blade is changed. This is shown in Figure 7-6. The effect of this change 
in angle of attack is to change the lift produced by the blade, which results in a non-uniform 
pressure wave - noise! . If the change in the angle of atack is large enough, the blade can 
stall and vibrate horribly. If the entire blade leading edge passes through this region of non- 
uniformity simultaneously, large unsteady pressure fluctuations result. In english, this means 
that a lot of noise is generated. If, however, the blade is skewed such that the blade sections 
individually experience the change in angle of attack, the strength of the unsteady pressure 
wave is reduced, and much less noise is generated. 

Reduced the pitch angle of the blade tip . The effect of the skew is to reduce the pitch angle 
of the tip section. This is important in fan applications where the meridional blade profile is a 
premium, and any reduction in the length of the meridional blade sections positively impacts 
the overall system design. 



The design process proceeded as follows: 

1. An initial-guess blade shape is designed with the desired skew distribution. In this case the 
distribution specified a tip skew of 40 degrees. 

2. By holding the leading edge line fixed, the propeller blade design code alters the position of 
the rest of the blade points to find the blade shape required to satisfy the kinematic boundary 
condition and produce the required loading. 
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Radial Circulation Distribution: Fan Design vs Analysis 




M.l.T. Marine Hydrodynamics Laboratory 



Figure 7-4: The design radial circulation distribution and final PBD14.4-MTFLOW analysis of the 
same blade. The blade was designed using PBD14.2 which uses an image representation of the hub 
and duct. Note the discrepancies in the tip and hub region. This is most likely due to the PBD14.4 
vortex lattice representation of the hub and duct. 



Ducted Fan Throughflow Domain 
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Figure 7-5: These two plots show the final MTFLOW inviscid solution streamlines. The plot on the 
right is a blow-up of the propeller tip region. Notice that the propeller leading edge tip operates in 
the extremely accelerated flow round the duct lip. This is a purely inviscid flow effect. 
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Design Velocity Triangle 
Angle of Attack 



Vapparent 



1 



Vin 



Angle of Attack as Blade 
Passes through Reduced 
Velocity Region 



Figure 7-6: As the blade passes through a region of inflow with reduced velocity, the velocity triangle 
clearly shows the accompanying increase in angle of attack. 

3. Because the meridional distribution of swirl has changed, this new swirl distribution is input 
into MTFLOW, and a new circumferentially mean (i.e., meridional) velocity distribution is 
found, which satisfies the fluid equations of motion. 

4. The updated velocity field is extracted from MTFLOW and input into the propeller blade 
design code, where a new blade shape is found. That is, the process loops back to Step 2. 

This process continues until the root mean square (rms) of the change in blade shape is less than 
a user-specified tolerance. 

In plain propeller designer English, the main effect of the blade design code is to find a suitable 
pitch distribution ( i.e set a local sectional angle of attack on the blade section), given the desired 
skew distribution. The camber distribution across the chord is specified by the user and can not be 
altered. 



7.5 40° Forward Skew Blade Design 

A forward skewed blade was designed in the same manner as the backward skewed blade. Starting 
from the desired circulation distribution, the blade was skewed such that the tip was skewed to a 
forward, or negative, angle of 40 degrees. 

7.6 Comparison Between the No Skew and Skewed Blades 

The plots shown in Figure 7-7 show the effects of skew upon the resultant blade shapes. The most 
noticable change is in the pitch angle of the blade tip region. 



Figure 7-8 dramatically highlights the pitch and camber redistributions introduced by the skewing 
of the blade. 

The only section properties not investigated were the chord and raked distribution. Unlike a 
marine propeller which needs long chord lengths to prevent cavitation, an airscrew uses longer chord 
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0 Degree Skew 





40 Degree Forward Skew 




Figure 7-7: The final fan designs for the back skew, no skew, and forward skew fans. The small plot, 
to the right of the fan shows the distribution of pitch (phi) and skew angle over the radius of the 
blade. 
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Distribution of Pitch Angle and Max Section Camber 
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Figure 7-8: The rake and max section camber for the final fan designs. Notice how the changes 
in section pitch ( i.e ., changes in local angle of attack) are accompanied by an offsetting change in 
section max camber. 



lengths to improve unsteady performance. The effects upon performance, specifically with regards 
to boundary layer separation, were not investigated in this thesis. 

Finally, rake is used by some designers as a complement to skew to introduce a non-uniform 
leading edge span line to any inflow unsteadiness. In the ducted fan case considered above, the 
application of rake to the de-pitched back skewed blade would enhance the fan performance, since 
the leading edge tip would be moved further away from the high velocity region of the duct lip. 



7.7 A Re- Analysis of the Skewed Blade Design 

A re-analysis of the skewed blade designs revealed that there may be problems with the PBD fan 
design mode. Specifically, the use of images to represent the duct in the design, led to an inaccurate 
flow representation in the tip region. With the highly skewed blades, the old paradigm of extending 
the image vortex lattice system radially outward from the tip breaks down, because a significant 
portion of the fan blade does not lie along a radial line exending inward from the tip. The end result 
is that the image vortex lattice system produces too weak an effect upon the rest of the blade, and 
the blade is overloaded. 

This problem is highlighted after a re-analysis of the designed for blade was done with a vortex 
lattice representation of the hub and duct. With a vortex lattice duct, the kinematic boundary 
condition is exactly enforced to produce the blade to blade variational flow. From Figure 7-9 we see 
that the blade designed using the image method is overloaded and will probably stall. 
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40 Degree Back Skew Fan Analysis with PBD 14.4 Vortex Lattice Hub and Duct 




M.l.T. Marin* Hydrodynamics Laboratory, 2/24/99 



Figure 7-9: The curve on the left represents the spanvvise loading distribution predicted by PBD 
using an image hub and duct representation. The curve on the right represents the spanvvise loading 
distribution predicted by PBD using a vortex lattice hub and duct representation. 
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7.7.1 The Back Skew on the Test Stand 



The 40 degree back skew fan was built and tested by Airflow Research. At the design speed, the 
fan stalled due to the extremely aggresive pitch of the blades. When the fan was powered down, a 
tremendous amount of airflow was produced, with an attendantly high torque. Although no data 
was actually collected, observation confirmed that the vortex lattice hub and duct predicted results 
were more in line with reality. 



7.8 Specific Problems of the PBD Fan Design Mode 

The largest, noted problem is the fact that PBD internally splines all flow quantities with respect to 
the axial location. In highly contracting flow domains, such as this ducted fan case, where there are 
nearly vertical boundaries, splming with respect to x leads to a multi-valued solution. In english. the 
resulting 3 rd order polynomial oscillates wildly. Any interpolations taken from such a polynomial 
will in no way represent the actual fluid properties at the interogated point in space. 

One solution to this problem is to spline with respect to the distance along a streamline. The 
difficulty then is how to interpolate to find the flow velocities at an arbitray point. For instance, for 
an arbitrary point, it is impossible to a prion know which streamline the point lies along, and the 
streamwise distance along the streamline at which the point lies. 

Another proposed solution is to transform the flow domain from an x — y coordinate system to a 
general curvi-linear £ — rj coordinate system. Then a spline with respect to £ will not be multi-valued. 

A second noted problem is in the pitch angle and number of vortex lattice elements placed 
upstream of the blade leading edge tip. When a the tip of the propeller has a low pitch angle, as in 
the case of the 40° back skew fan, the current implementation of making the vortex lattice follow 
the pitch line of the blade leading edge results in a highly skewed lattice. 
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Chapter 8 

Conclusions 



MTFLOW coupled with PBD works. This thesis has shown, through the nominal and effective 
wake calculation comparison, that the coupling is based upon sound principles, and that there are 
no fundamental reasons that PBD and MTFLOW should not be coupled. As such, MTFLOW 
represents a viable alternative to axisvmmetric RANS codes for solving circumferentially mean 
throughflow problems, because of the large reduction in computational cost. 

The devil, though, is in the details. More refinement is needed in handling the streamline radial 
redistribution problem. But now that the coupling codes exist, these are refinements upon an 
established base. 



8.1 Propeller Blade Design Code Improvements 

The steady design and analysis propeller code, PBD, should be continually updated to reflect the 
increasing sophistication of preliminary designs. Because the blade singularity discretization leads to 
a requirement to couple PBD with some sort of throughflow solver, time should be spent on creating 
common data models which can easily be manipulated by PBD and any generic throughflow solver 
such that information is not lost in the translation. 

8.1.1 Analysis Mode 

The analysis capabilities of PBD are extremely robust. The extended vortex lattice hub and duct 
system, however, increases the computational time associated with the solution. This increase in 
time becomes noticable when coupling PBD with an extremely fast throughflow solver such as 
MTFLOW. More attention to the underlying numerical data representation, and numerical solution 
techniques, with an eye towards speed, would propel PBD back ahead of the pack. 

There may be issues with the current implementation of the hub and duct vortex lattice model. 
One issue is the increase in computational time required due to the larger matrix system involved. 
Another issue is the lack of user control over the orientation and placement of the lattice grid. Highly 
stretch grid lattices may result from various geometries. Perhaps a solution to these problems may 
be a B spline hub and duct representation. 

8.1.2 Design Mode 

Attention should be focused on incorporating a vortex lattice hub and duct in the design mode to 
account for the effects of blade skew and radially non-uniform hub and ducts. Perhaps the only 
solution will become, then, design by analysis. Perhaps, though, there is the elegant solution lurking 
just around the corner. 
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8.2 Boundary Layer Modelling Improvements 

This thesis highlighted the disimilar boundary layer profiles which result from the use of present 
state of the art RANS and IBLT boundary layer models. The fact is that the question is still open on 
if a generally applicable, turbulent boundary layer model exists. An added dimension of complexity 
is the proper treatment of the propeller induced body forces in the turbulent boundary layer closure 
equations. Such research may lead into a better understanding of flow around submerged bodies, 
and better predict the propeller inflow velocity field. 
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Appendix A 



User Notes on Running MTFLOW 



A.l Problem Solving 

When problems are encountered wit h MTFLOW, it is best to look at the contour plot within MTSOL 
of the state parameters /?, </, and P. It is also useful to look the MTSOL contour plot of any variables 
that you input to check for spline problems. Spline problems manifest themselves as large or small 
holes in the solution. I found that the velocity cont.out plot (MTFLOW q variable) was helpful to 
identify trouble areas. 



A. 2 Grid Issues 

The interpolation of the input swirl into MTFLOW is somewhat dependent upon the grid. While 
arbitrarily increasing grid density will not usually solve interpolation problems, clustering t lines 
near the problem area can sometimes fix the problem. 



A. 3 Input Swirl: rAV# 

A. 3.1 Interpolation 

It is critical to adequately refine the swirl input near the blade tip. This is increasingly important 
when several radii are redistributed by PBD2MT near the tip. This require bunching several lines 
of constant t near the tip radius. 

A.3.2 Convection 

The user must manually convect any excess swirl downstream. If this downstream convected swirl 
encounters a downstream mechanical device that adds or removes swirl, the operator must combine 
the upstream and local swirl in a logical manner. The solution that I used was to first build the 
tflow input file for the downstream blade row in the usual manner and then convect the swirl from 
the trailing edge of the upstream blade row along the previous iteration streamlines through the 
downstream blade row. Remember that as the swirl convects down streamlines that change radial 
distance off the centerline of the machine, the r&Vg must be changed accordingly, since internal to 
MTFLOW the V t heta is calculated by dividing this quantity by the local radius. 
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A. 4 Open Propellers 

A. 4.1 Far Field Setting 

Assume a constant pressure jet boundary and set the Far Field singularity equal to 0. Inside 
MTSOL, set F F — 3. While this may cause a slight degradation of the flow boundary condition, 
my experience has shown that using a singularity causes more problems, especially when modelling 
a domain with an upstream or downstream shaft. 

A. 5 The Mach Number in Water 

When running MTFLOW, the user specifies a flow Mach number and Reynold’s number. While the 
Reynold’s number is widely used in marine fluid dynamics, the use of mach number, in a nominally 
incompressible fluid, warrants a moment of thought. For a general case of interest, start from the 
known Reynold’s number of the problem. 




v 



( A. 1 ) 



Knowing the length scale L and assuming a kinematic viscosity, the velocity scale is determined. 



V = 



Rev 

~T 



(A. 2) 



If it is now assumed that the speed of sound in water is nominally 1680 ™ [25], the mach number 
for the flow problem is known. 



M = 



Rev 
1680 L 



(A.3) 



Since the speed of sound in MTFLOW, designated as a in MTFLOW, is always 1.0, the mach 
number calculated from equation (A.3) is directly input into MTFLOW as the inlet mach number. 
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Appendix B 



Modifications to MTFLOW Source 
Code 



Modifications to the MTSET/MTSOL/MTFLO family of programs was limited to the mtset.f and 
io.f subroutines. 



B.l Modifications to mtset.f 

Within mtset.f. an output routine was added to output the original gridlines calculated by MTSET. 
This information is required when BL2BODY lays the calculated boundary layer velocity profile 
back into the flow domain. 



B.2 Modifications to io.f 

The io.f subroutine was modified to output the integral boundary layer quantities and the flow 
solution with grid streamline positions. It has been validated for both external propellers and 
internal flow cases. 

The following output files are written. 



SUBROUTINE 


ALTERATIONS 


mtset.f 


output ORIGGRID.tec, original grid set by MT- 
SET, over the solid wall body points (downstream 
of the leading edge) 


io.f 


output OUTGRID.tec, final grid from MTSOL 
solution 

output OUTVEL.tec, cell-centered velocities from 
MTSOL 

output OUTBL.tec, boundary layer quantities in 
boundary layer coordinates 



Table B.l: Custom MTFLOW output file in tecplot formatted ascii for tecplot viewing and use by 
BL2BODY 
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SUBROUTINE 


ALTERATIONS 


io.f 


output OUTGRID.tec, final grid from MTSOL so- 
lution. It is useful for ensuring that BL2B0DY 
routines did a good job in accurately reconstruct- 
ing to flow domain. 

output OUTVEL.tec, This file has the c, r loca- 
tion for the Ai'V# velocity components. It is 

cell- centered velocity data, like most RANS finite 
volume solver output data. 

output OUTBL.tec, boundary layer quantities in 
boundary layer coordinates 



Table B.2: Custom MTFLOW output file in tecplot formatted ascii for tecplot viewing and use by 
BL2B0DY 
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Appendix C 



ITTC Propeller 4119 Input Files 



C.l PBD Input File 

pbd-14.4 analysis of p4119 - coupled 

P4119.BSN 

restart . vel 

3 15 15 
1 

14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 

7 0.0 0 0.0 

8 6 

0.84 300 
0.3 0.3 
10 -10 0 0 
6 
0 

1 0.001 0.00 0.0 1 
1 0.02 

4 2 

1.11023 1.000 1.500 0.02 

0.00000 0.0000 0.0000 0.0000 0.0000 
0.20000 0.3000 0.4000 0.5000 0.6000 
0.06576 0.0563 0.0477 0.0396 0.0321 
0.00814 0.0086 0.0085 0.0084 0.0082 
0.000 0.000 0.000 0.000 0.000 0.000 
0.000 0.000 0.000 0.000 0.000 0.000 
0.000 0.000 0.000 0.000 0.000 0.000 
0.000 0.000 0.000 0.000 0.000 0.000 



MTFLOW analysis 

: FILE NAME FOR BLADE B-SPLINE NET 
: FILE NAME FOR WAKE FIELD 
: nblade nkey , mkey 
: ispn (0=unif orm, l=cos) 

15 16 17 18 19 :mctrp ,mc(n) 

: ihub,hgap, iduc ,dgap 
: HDWAK , NTWAKE 
: Cq, MXITER 
: OVHANG(l), 0VHANG(3) 

: nx ,ngcoef f ,mltype ,mthick 
: imode 
: nwimax 

: niter , tweak, bulge ,radwgt ,nuf ix 



: nplot,hubshk 
: NOPT, NBLK 

: ADVCO XULT XFINAL DTPROP 
0.0000 0.0000 0.0000 0.0000 0.000 : G 

0.7000 0.8000 0.9000 0.9500 1.0000 : r/R 

0.0250 0.0183 0.0120 0.0089 0.00000 : T/D 

0.0081 0.0080 0.0079 0.0078 0.0079 : CD 

0.000 0.000 0.000 0.000 : ua 

0.000 0.000 0.000 0.000 : uau 

0.000 0.000 0.000 0.000 : ut 

0.000 0.000 o.ooa 0.000 : utu 



C.2 MTFLOW Input Files 


C.2.1 Coupling Admin File 


432000 


Reynolds number 


0.001488 


inlet mach number 


1.0093 


Vship used in PBD to find J 


15 


aft most point of body 


4.000 


x location of LE tip 


1.000 


r location of LE tip 


tf low .4119 


tflow file name 



75 



walls .4119 
p41 19 . pbd 



! walls file name 
! pbd input file name 



C.2.2 Walls. 4119 File 

4119 ITTC walls case 

1.000000 7.500000 0.000000 2.000000 



7.600000 


0.200000 


5.449084 

5.394665 

5.340678 

5.287190 


0.200000 

0.200000 

0.200000 

0.200000 


5.234253 

5.181955 

5.130387 

5.079632 

5.029772 


0.200000 

0.200000 

0.200000 

0.200000 

0.200000 


4.980885 


0.200000 


4.933049 


0.200000 


4.886336 

4.840815 

4.796552 


0.200000 

0.200000 

0.200000 


4.753613 

4.712055 


0.200000 

0.200000 


4.671937 


0.200000 


4.633312 


0.200000 


4.596230 


0.200000 


4.560740 

4.526885 


0.200000 

0.200000 


4.494707 

4.464243 


0.200000 

0.200000 


4.435529 


0.200000 


4.408596 


0.200000 


4.383472 


0.200000 


4.360183 

4.338752 


0.200000 

0.200000 


4.319196 


0.200000 


4.301532 

4.285772 

4.271927 


0.200000 

0.200000 

0.200000 


4.260002 


0.200000 


4.250000 


0.200000 


4.240000 


0.200000 


4.230000 


0.200000 


4.220000 


0.200000 


4.210000 

4.200000 


0.200000 

0.200000 


4.190000 


0.200000 


4.180000 

4.170000 


0.200000 

0.200000 


4.160000 


0.200000 


4.150000 

4.140000 


0.200000 

0.200000 


4.130000 


0.200000 
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4.120000 

4.110000 

4.100000 

4.090000 

4.080000 

4.070000 

4.060000 

4.050000 

4.040000 

4.030000 

4.020000 

4.010000 

4.000000 

3.250000 

3.237500 

3.225000 

3.212500 

3.200000 

3.188016 

3.176374 

3.165110 

3.154251 

3.143820 

3.133834 

3.124304 

3.115236 

3.106635 

3.098497 

3.090819 

3.083593 

3.076809 

3.070454 

3.064515 

3.058977 

3.053822 

3.049034 

3.044594 

3.040485 

3.036688 

3.033185 

3.029959 

3.026990 

3.024263 

3.021761 

3.019467 

3.017367 

3.015447 

3.013693 

3.012093 

3.010637 

3.009314 

3.008113 

3.007026 

3.006046 



0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.199641 

0.198600 

0.196933 

0.194697 

0.191947 

0.188738 

0.185122 

0.181149 

0.176870 

0.172329 

0.167569 

0.162633 

0.157556 

0.152374 

0.147119 

0.141818 

0.136499 

0.131184 

0.125893 

0.120644 

0.115453 

0.110331 

0.105290 

0.100338 

0.095481 

0.090723 

0.086069 

0.081518 

0.077072 

0.072730 

0.068492 

0.064357 

0.060322 

0.056386 

0.052547 

0.048802 



3.005163 

3.004371 

3.003665 

3.003037 

3.002482 

3.001996 

3.001573 

3.001210 

3.000901 

3.000645 

3.000436 

3.000271 

3.000148 

3.000064 

3.000016 

3.000000 

3.000016 

3.000064 

3.000148 

3.000271 

3.000436 

3.000645 

3.000901 

3.001210 

3.001573 

3.001996 

3.002482 

3.003037 

3.003665 

3.004371 

3.005163 

3.006046 

3.007026 

3.008113 

3.009314 

3.010637 

3.012093 

3.013693 

3.015447 

3.017367 

3.019467 

3.021761 

3.024263 

3.026990 

3.029959 

3.033185 

3.036688 

3.040485 

3.044594 

3.049034 

3.053822 

3.058977 

3.064515 

3.070454 



0.045150 

0.041587 

0.038112 

0.034722 

0.031414 

0.028186 

0.025038 

0.021965 

0.018968 

0.016045 

0.013194 

0.010414 

0.007705 

0.005067 

0.002498 

0.000000 

- 0.002498 

- 0.005067 

- 0.007705 

- 0.010414 

- 0.013194 

- 0.016045 

- 0.018968 

- 0.021965 

- 0.025038 

- 0.028186 

- 0.031414 

- 0.034722 

- 0.038112 

- 0.041587 

- 0.045150 

- 0.048802 

- 0.052547 

- 0.056386 

- 0.060322 

- 0.064357 

- 0.068492 

- 0.072730 

- 0.077072 

- 0.081518 

- 0.086069 

- 0.090723 

- 0.095481 

- 0.100338 

- 0.105290 

- 0.110331 

- 0.115453 

- 0.120644 

- 0.125893 

- 0.131184 

- 0.136499 

- 0.141818 

- 0.147119 

- 0.152374 



3.076809 

3.083593 

3.090819 

3.098497 

3.106635 

3.115236 

3.124304 

3.133834 

3.143820 

3.154251 

3.165110 

3.176374 

3.188016 

3.200000 

3.212500 

3.225000 

3.237500 

3.250000 



- 0.157556 

- 0.162633 

- 0.167569 

- 0.172329 

- 0.176870 

- 0.181149 

- 0.185122 

- 0.188738 

- 0.191947 

- 0.194697 

- 0.196933 

- 0.198600 

- 0.199641 

- 0.200000 

- 0.200000 

- 0.200000 

- 0.200000 

- 0.200000 
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Appendix D 

Extreme Computations 

This chapter represents output from PBD and MTFLOW that may indicate problems with the 
solution. The solutions are also included to guide the curious reader. 



D.l MTFLOW Spline Problems 

Because the input data is splined within MTFLOW, any non-smooth data results in spline interpo- 
lation problems. If these problems are ignored, unrealistic flow fields like the one shown in Figure 
D-l result. 



D.2 PBD Wake Routine Problems 

A viscous boundary layer along the downstream hub or casing (duct) can lead to problems. Because 
the meridional fluid velocity is slower at the hub and tip radii and faster in the propeller jet, the 
wake will appear to bow in the middle. 



The solution is to improve the velocity smearing scheme used within the program VELCON such 
that the meridional velocity field is more nearly uniform. 
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Axial Velocity Contours from MTFLOW 



Z 




Figure D-l: The flowfield resulting from an MTFLOW input with spline interpolation problems. 
This physically unrealistic flowfield is handed back to PBD, which solves the blade analysis problem. 
The results are highly suspect, even though convergence may be achieved. 
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Figure D-2: Wake shot of 3-bladed 4119 with viscous boundary layer on the shaft. 
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